Construction and integration of genetic linkage maps from three multi-parent advanced generation inter-cross populations in rice

Background The construction of genetic maps based on molecular markers is a crucial step in rice genetic and genomic studies. Pure lines derived from multiple parents provide more abundant genetic variation than those from bi-parent populations. Two four-parent pure-line populations (4PL1 and 4PL2) and one eight-parent pure-line population (8PL) were developed from eight homozygous indica varieties of rice by the International Rice Research Institute (IRRI). To the best of our knowledge, there have been no reports on linkage map construction and their integration in multi-parent populations of rice. Results We constructed linkage maps for the three multi-parent populations and conducted quantitative trait locus (QTL) mapping for heading date (HD) and plant height (PH) based on the three maps by inclusive composite interval mapping (ICIM). An integrated map was built from the three individual maps and used for QTL projection and meta-analysis. QTL mapping of the three populations was also conducted based on the integrated map, and the mapping results were compared with those from meta-analysis. The three linkage maps developed for 8PL, 4PL1 and 4PL2 had 5905, 4354 and 5464 bins and were 1290.16, 1720.01 and 1560.30 cM in length, respectively. The integrated map was 3022.08 cM in length and contained 10,033 bins. Based on the three linkage maps, 3, 7 and 9 QTLs were detected for HD while 6, 9 and 10 QTLs were detected for PH in 8PL, 4PL1 and 4PL2, respectively. In contrast, 19 and 25 QTLs were identified for HD and PH by meta-analysis using the integrated map, respectively. Based on the integrated map, 5, 9, and 10 QTLs were detected for HD while 3, 10, and 12 QTLs were detected for PH in 8PL, 4PL1 and 4PL2, respectively. Eleven of these 49 QTLs coincided with those from the meta-analysis. Conclusions In this study, we reported the first rice linkage map constructed from one eight-parent recombinant inbred line (RIL) population and the first integrated map from three multi-parent populations, which provide essential information for QTL linkage mapping, meta-analysis, and map-based cloning in rice genetics and breeding.


Background
Rice is an important staple-food crop for nearly half of the world's population. According to the International Grain Council (IGC), global rice production is projected to reach 505 million tons in 2019-2020, approximately 90% of which will be produced and consumed by Asians (Kong et al. 2015;Bazrkar-Khatibani et al. 2019). In the last several decades, the rice yield has increased slowly in most production countries, and breaking the barrier to yield growth has become a difficult challenge for breeders (Peng et al. 2004;Meng et al. 2016). It is promising to break through the yield stagnation by integrating conventional breeding methods with molecular marker technology and genomics (Liang et al. 2016).
Construction of linkage map by using molecular markers is a crucial step in genetic and genomic analysis that facilitates the discovery of quantitative trait locus/ loci (QTLs) and provides chromosomal information for gene cloning and marker-assisted breeding (Meng et al. 2015;Palumbo et al. 2019). Since the publication of the first rice linkage map, which included 135 restriction fragment length polymorphism (RFLP) markers (McCouch et al. 1988), great progresses have been made in rice linkage map construction, especially with the development of high-throughput sequencing technology in recent years. Researchers have constructed a large quantity of maps in rice by using different types of populations and various molecular markers. For example, Harushima et al. (1998) reported a genetic linkage map containing 2275 cDNA markers covering 1521.6 cM in an F 2 population of Nipponbare × Kasalath. Yin et al. (2015) constructed a linkage map containing 143 SSR markers in one japonica × indica genetic population consisting of 215 recombinant inbred lines (RILs). De Leon et al. (2016) constructed a high-density GBSderived SNP linkage map with 2817 bins from a rice RIL population. Swamy et al. (2018) published 6 K SNP chipderived SNP linkage maps in two doubled-haploid (DH) populations of rice. Kim (2018) developed a SNP map with 1954 bins for an F 2 population.
Low marker number and density, and large number of co-localized markers are limitations of linkage maps constructed from single mapping populations, whereas consensus maps combining genetic information from multiple populations provide better genome coverage with higher marker density (Galeano et al. 2011). A consensus map can also validate marker order, identify chromosomal rearrangements and confirm the position of common markers across mapping populations (Salvi et al. 2009;Maccaferri et al. 2015). Furthermore, comparisons of QTLs on different linkage maps are based on the consensus map (Wen et al. 2017). Some tools for consensus map construction have been developed and widely used, such as JoinMap (Van Ooijen 2006), MergeMap (Wu et al. 2010), BioMercator (Arcade et al. 2004), LPmerge (Endelman and Plomion, 2014) and MultiPoint (Ronin et al. 2012). Using these tools, some integrated maps from bi-parent populations have been reported in rice.  integrated 15 individual maps, resulting in a consensus map with 531 markers and a total map length of 1821 cM. Swamy and Sarla (2011) constructed a consensus map from 11 maps with 699 markers and a map length of 1676 cM. Wu et al. (2016) developed a consensus map from 6 maps with 6970 markers and a total length of 1823.1 cM. Lei et al. (2018) constructed consensus maps from 6 maps with 7156 markers and a length of 1112.71 cM. Islam et al. (2019) developed a high-density consensus map with 12,327 simple sequence repeat (SSR) and SNP markers and a length of 2936 cM.
In the past 15 years, multi-parent advanced generation inter-cross (MAGIC) populations have become increasingly popular and have been developed in various plant species (Cavanagh et al. 2008;Kover et al. 2009;Huang et al. 2012;Mackay et al. 2014;Dell'Acqua et al. 2015;Pascual et al. 2015;Sannemann et al. 2015). Some MAGIC populations have been reported in rice. For example, International Rice Research Institute (IRRI) developed four MAGIC populations in rice (Bandillo et al. 2013), andOgawa et al. (2018) developed a MAGIC population from eight high-yielding rice cultivars, and GWAS-based QTL mapping studies have also been conducted in those populations (Shen et al. 2017;Ogawa et al. 2018). Compared with bi-parent populations, MAGIC populations contain multiple alleles and possess an amplified number of recombination events and greater genetic diversity. However, the larger numbers of alleles and marker types at each locus in MAGIC populations complicate recombination frequency estimation and linkage map construction. For most MAGIC populations, genome-wide association studies (GWAS) have been employed for gene detection without linkage maps. However, GWAS have some disadvantages, such as lack of background control, and it doesn't take advantage of linkage information. Accurate linkage analysis and linkage QTL mapping are essential for MAGIC populations. Several R packages and software have been developed to provide functions for constructing linkage maps and QTL mapping in MAGIC populations, including R/qtl (Broman et al. 2003), R/happy (Mott et al. 2000), R/mpMap (Huang and George 2011), and GAPL (integrated genetic analysis software for multi-parent pure-line populations) (Zhang et al. 2017;Zhang et al. 2019). Among these packages, GAPL is a stand-alone package used to build accurate linkage maps in a short time for MAGIC populations . Inclusive composite interval mapping (ICIM) implemented in GAPL is applicable for QTL mapping in both four-parent and eight-parent pure-line populations, with high detection power, a low false discovery rate and accurate estimation of QTL positions and effects (Shi et al. 2019).
To further facilitate QTL discovery and variety development in rice, IRRI developed three rice MAGIC populations from eight diverse indica rice varieties (Meng et al. 2016), i.e., two four-parent RIL populations and one eight-parent RIL population. Meng et al. (2016Meng et al. ( , 2017 and Shen et al. (2017) conducted GWAS of abiotic tolerance, heading date (HD) and plant height (PH) in these populations. In this study, we focused on construction of linkage maps and integrated map, and linkage QTL mapping for the three multi-parent populations. Our objectives were (1) to construct individual linkage maps for the three multi-parent populations; (2) to construct an integrated map from the three individual maps; and (3) to conduct QTL mapping of HD and PH using the integrated map and compare the mapping results with those from individual populations with metaanalysis.
Field trials were conducted at the Pingxiang Experimental Station (113.85°N, 27.6°E) in Jiangxi Province (Shen et al. 2017). Seeds were sown in a seedling nursery on May 30, 2016. Seedlings were transplanted into the field 25 days after sowing. The three RIL populations, including 392 lines of 4PL1, 386 lines of 4PL2 and 566 lines of 8PL were planted in a randomized complete block design with two replications, with 10 plants in a row for each replication. The plant spacing was 17 cm between rows and 20 cm within each row. HD was recorded as the number of days from sowing to the date when 50% of the plants headed, whereas PH was measured based on five representative plants after complete heading. Mean values from the two replications were used for statistical analysis and QTL mapping.

DNA Extraction and Genotyping
Lines from the three MAGIC populations were first genotyped and reported by Meng et al. (2017). The eight parents were re-genotyped in this study to obtain more accurate genotypic data and reduce the rate of missing data. Genomic DNA for SNP genotyping was isolated from approximately 100-mg fresh leaf samples of 5-weekold seedlings using a modified cetyltrimethylammonium bromide (CTAB) method (Murray and Thompson, 1980). Genotyping was performed using a customized rice 55 K SNP array containing 56,897 SNP screened from the 3 K Rice Genome Project (Project 2014; Meng et al. 2017;Wang et al. 2018). Target DNA preparation, chip hybridization, and array processing were conducted by CapitalBio Technology (Beijing, China) according to the Affymetrix Axiom® 2.0 assay protocol. Finally, 39,070 high-quality (poly high-resolution) SNPs with genetic diversity were selected based on comparison to the Oryza sativa cv. Nipponbare IRGSP 1.0 reference genome (Kawahara et al. 2013). For convenience, SNPs were renamed based on their physical locations.

Construction of Individual and Integrated Maps
Quality control of the genotypic data was performed for the three populations separately using GAPL V1.2 software . First, using the SNP (SNP genotypic data conversion) functionality in GAPL, DNA base data were converted into a format that could be imported into the software. SNPs showing nonpolymorphism in parents or progenies or that were missing in one or more parents were filtered out. Second, using the BIN (binning of redundant markers) function, markers with more than 10% missing values were deleted, meanwhile, the marker with the minimum percentage of missing data was selected to represent the co-localized markers. A set of co-localized markers in a specific population was defined as one bin. Finally, markers having more than 3% heterozygosity in progenies and progenies having more than 5% heterozygosity were removed. SNPs retained after filtering were used for map construction.
Linkage maps of 8PL, 4PL1 and 4PL2 were constructed separately using the PLM (map construction in multi-parent derived pure-line populations) function in GAPL. First, SNPs with a known chromosome ID on a physical map were anchored. The unanchored SNP markers were grouped based on a minimum logarithm of odds (LOD) score of 3. Second, the nearest-neighbour algorithm and two-opt algorithm from the travelling salesman problem (TSP, Lin and Kernighan, 1973) were adopted for marker ordering. Finally, the sum of adjacent recombination frequencies calculated with a five-marker window size was used as the rippling criterion. Recombination frequency was converted into map distance by the Kosambi mapping function.
The integrated map was built by using MergeMap V4.1 software (Wu et al. 2010) based on the individual maps from 8PL, 4PL1 and 4PL2. Equal weight was given to the three individual maps. The R package Linkage-MapView was used for visualization of the constructed maps (Ouellette et al. 2017).

QTL Analysis of the Three Populations
QTL mapping of HD and PH was conducted separately in 8PL, 4PL1 and 4PL2 by using ICIM with the PLQ (QTL or gene detection in multi-parent derived pureline populations) function in GAPL. The scanning step was set at 0.1 cM. Probabilities of adding and removing variables in stepwise regression were set at 0.001 and 0.002, respectively. The LOD threshold was determined by a permutation test with 1000 runs, and the type-I error was set to 0.05. Mapping results from 8PL were compared with the results obtained by GWAS in Shen et al. (2017), where 0.00001 was used as the threshold for the P-value.

QTL Projection and Meta-QTL Analysis
BioMercator V4.2 was used to project the detected HD and PH QTLs onto the integrated map and perform QTL meta-analysis (Veyrieras et al. 2007, Sosnowski et al. 2012. Based on the estimated position, LOD score, phenotypic variation explained (PVE) by each QTL, and confidence interval, these QTLs were projected onto the integrated map. Meta-analysis of the projected QTL clusters was performed using Goffinet's and Gerber's algorithm (Goffinet and Gerber 2000) for each chromosome. The lowest Akaike information criterion (AIC) value was used to select the best QTL model. The position and 95% confidence interval of each meta-QTL (MQTL) were calculated. QTL mapping of DH and PH using the PLQ function in GAPL was also conducted with the integrated map. The mapping parameters were the same as those used for QTL mapping with individual maps.

General Information on the Genotypic and Phenotypic Data
A total of 39,070 markers were used to filter the pipeline in the three MAGIC populations. The missing rate of these markers is shown in Additional file 8: Figure 1A.
In 8PL, 4PL1 and 4PL2, there were 279, 252 and 270 markers whose missing rates were greater than 10%; percent of 13.17, 12.03, and 4.48 of the markers had a heterozygosity level over 3%, respectively (Additional file 8: Figure 1B). After removing these markers as well as SNPs with non-polymorphism in parents or progenies or that were missing in one or more parents, 13,254, 10, 825 and 13,479 SNPs were retained for 8PL, 4PL1 and 4PL2, respectively. A total of 518, 364 and 356 progenies were kept for 8PL, 4PL1 and 4PL2, respectively, after controlling for progeny heterozygosity. After removing co-localized markers, 5906, 4368 and 5481 markers were retained in the 8PL, 4PL1 and 4PL2 populations and then used for linkage map construction, respectively. Populations of 8PL and 4PL1 had 2352 common markers; populations of 8PL and 4PL2 had 2147 common markers; populations of 4PL1 and 4PL2 had 1379 common markers; and the three populations had 942 markers in common (Additional file 8: Figure 1C). The phenotypic variation in the eight parents and the three RILs populations are shown in Additional files 1 and 2: Tables S1 and S2, respectively. The range of variations in progenies was greater than that in parents, reflecting the presence of transgressive segregation. The kurtosis and skewness of HD and PH were close to zero, except for those of PH in 4PL2. The broad-sense hereditability (H 2 ) was 0.9867, 0.9346 and 0.9736 for HD and 0.7766, 0.7711 and 0.9217 for PH in 8PL, 4PL1 and 4PL2, respectively.

Individual and Integrated Maps
The distribution of SNPs and linkage map information for the three populations are shown in Table 1. For 8PL, the full genome was 1290.16 cM in length and included 5906 markers across the 12 chromosomes (Fig. 2). The number of bins was 5905, i.e., 5905 markers had unique map positions. The average distance between adjacent bins was 0.22 cM. The most saturated chromosome was chromosome 8, which had an average bin distance of 0.15 cM, while chromosome 6 had the largest average bin distance of 0.39 cM. The average chromosome length was 107.51 cM, with the longest being chromosome 1 (150.66 cM) and the shortest being chromosome 12 (67.79 cM) ( Table 1). For 4PL1, the entire genome spanned 1720.01 cM and contained 4368 markers (Fig. 3). The number of bins was 4354. The average bin distance was 0.40 cM, and the average chromosome length was 143.35 cM. Chromosome 7 had the smallest average bin distance of 0.28 cM, whereas chromosome 6 had the largest average bin distance of 0.66 cM. Chromosome 1 was the longest, i.e., 201.39 cM in length, and chromosome 10 was the shortest, i.e., 100.54 cM in length. For 4PL2, the entire genome was 1560.30 cM in length and contained 5481 markers (Fig. 4). The number of bins was 5464. The average bin distance was 0.29 cM, with the smallest average bin distance of 0.22 cM on chromosome 4, whereas the largest average bin distance of 0.44 cM on chromosome 3. The average chromosome length was 130.03 cM, with the longest being chromosome 1 (207.53 cM) and the shortest being chromosome 10 (82.39 cM). Markers were approximately evenly distributed among the three individual maps from the three populations. The longest bin intervals were 5.64, 12.98, and 13.14 cM for 8PL, 4PL1 and 4PL2, respectively, whereas the shortest bin intervals were all equal to 0.1 cM. The 8PL map had the largest marker number, shortest genome length, and shortest average bin distance, whereas the 4PL1 map had the smallest marker number, longest genome length, and longest average bin distance (Table 1, Figs. 2, 3 and 4).
Contribution of each parent can be quantified by the proportion of the parental genome in the progeny (Wang and Bernardo 2000). In MAGIC population, it cannot be directly estimated by the bi-allelic SNP markers. However, after the map construction in GAPL, any missing and incomplete marker (e.g. bi-allelic SNP marker) can be imputed as complete marker, i.e., locus with eight alleles. Therefore, contribution of each parent can be calculated after imputation. In 4PL1, the proportion was 0.3183, 0.1184, 0.2683 and 0.2950 for parents A to D, respectively. In 4PL2, the proportion was 0.3536, 0.1725, 0.2318 and 0.2421 for parents E to H, respectively. In 8PL, the proportion was 0.1372, 0.1043, 0.1314, 0.1198, 0.1901, 0.0696, 0.1236 and 0.1240 for parents A to H, respectively. Theoretically, each parent contributed equally to MAGIC populations. The variation on parental contribution observed may be caused by genetic drift, potential selection, marker ascertainment bias and sampling errors.
The integrated map based on the three individual maps comprised 10,785 markers, 10,033 of which had unique map positions (Fig. 5). The genome was 3022 cM in length, with an average distance of 0.30 cM between adjacent bins (Table 2). Chromosome 1 had the largest number of 1403 SNPs, and chromosome 10 contained the smallest number of 639 SNPs. The average chromosome length was 251.84 cM. Chromosome 1 was the longest chromosome with total length of 385.18 cM and a mean distance of 0.29 cM between adjacent bins. Chromosome 10 was the shortest with total length of 139.54 cM and a mean bin distance of 0.24 cM. Chromosome 10 had the smallest average bin size of 0.24 cM, whereas chromosome 6 had the largest average bin size of 0.41 cM. The integrated map had no gaps longer than 17.05 cM. The average bin distance in the integrated map was larger than that in 8PL but smaller than that in 4PL1 and 4PL2. The genome length of the integrated map was greater than the three individual maps. Colocalized markers in some individual maps could be separated in other maps. Those markers will also be separated in the integrated map, leading to a longer genome in integrated map.

QTL Analysis Based on the Three Linkage Maps
From permutation tests, LOD thresholds at the 0.05 significance level were calculated at 8.59, 5.22 and 5.40 for 8PL, 4PL1 and 4PL2, respectively. The estimated positions and effects of the detected QTLs are given in Table 3 for 8PL and Additional files 3 and 4: Tables S3 and S4 for 4PL1 and 4PL2, respectively.
For HD, three QTLs were detected in 8PL, explaining a total of 24.29% of the phenotypic variance. Two were located on chromosome 6, and one, on chromosome 8. qHD6.1, located at 70.50 cM on chromosome 6, had the largest LOD score of 42.38 and the largest PVE of 15.66% (Table 3). At qHD6.1, parents A, B, C, E, and F provided alleles that reduced HD. Given that the length of the confidence interval was 0.5 Mb, i.e., the distance between a gene and an adjacent marker of the detected QTL was less than 0.5 Mb on the rice physical map, the gene and the QTL were declared to coincide. qHD6.2 and qHD8 coincided with Hd3a (LOC_Os06g06320), located from 2,940,004 to 2,942,452 bp (Taoka et al. 2011), and Ghd8 (LOC_Os08g07740), from 4,334,739 to 4,333,846 bp (Yan et al. 2011), respectively. In 4PL1, seven HD QTLs were detected, explaining a total of 60.18% of the phenotypic variance. Two were located on chromosome 8, and one each, on chromosomes 3, 4, 6, 7 and 11 (Additional file 3: Table S3). qHD6, located at 60.9 cM on chromosome 6 and coinciding with Hd3a, had the largest LOD score of 51.86 and the largest PVE of 28.63%. In 4PL2, nine QTLs affecting HD were detected, explaining a total of 55.83% of the phenotypic variance, among which qHD8 had the largest LOD score of 54.71 and largest PVE of 25.07% and was close to Ghd8 (Additional file 4: Table S4). In addition, qHD3.1, located at 13.10 cM on chromosome 3, was close to the reported gene HD6 (LOC_Os03g55389) that was located from 31,  (Yamamoto et al. 2000). For PH, six QTLs were detected in 8PL, explaining a total of 19.68% of the phenotypic variance. Two were located on chromosome 1, and one each, on chromosomes 2, 8, 11 and 12 (Table 3). qPH1.1 was located at 7.10 cM on chromosome 1, having the largest LOD score of 44.81 and the largest PVE of 6.92%. Alleles from parents A, B, D, E, G and H at qPH1.1 reduced PH. The cloned dwarf gene sd1 (LOC_Os01g66100), located from 38,381,423 to 38,384,165 bp (Monna et al. 2002), and the semi-dwarf gene Psd1 (LOC_Os01g60740), located from 35,129,858 to 35,130,917 bp (Li et al. 2014), were close to qPH1.1 and qPH1.2. Nine QTLs were detected in 4PL1, two each on chromosomes 5 and 8 and one each on chromosomes 1, 3, 6, 7 and 10, explaining a total of 74.99% of the phenotypic variance (Additional file 3: Table S3). qPH1 was located at 181.7 cM on chromosome 1, near the cloned gene sd1, and had the largest LOD score of 98.95 and largest PVE of 55.80%. In 4PL2, 10 QTLs were detected, located on chromosomes 1, 2, 5, 7, 8, and 10 and explaining a total of 42.33% of the phenotypic variance (Additional file 4: Table S4). qPH1.2, with the largest LOD score of 164.76 and the largest PVE of 34.85, coincided with d61, located from 29,927,543 to 29,931,487 bp (Yamamuro et al. 2000).

QTL Detection Based on Meta-Analysis
Based on the QTLs detected in 8PL, 4PL1 and 4PL2, a total of 44 MQTLs were identified on the integrated map, including 19 for HD and 25 for PH (Table 4, Additional file 9: Figure S2). Four of them were detected in two populations, i.e. MqHD3.2, MqHD11.2 and MqPH10.2 was detected in both 4PL1 and 4PL2; MqPH8.2 was detected

QTL Analysis Based on the Integrated Map
From permutation tests, LOD thresholds at the 0.05 significance level were calculated at 9.08, 5.71 and 5.87 for 8PL, 4PL1 and 4PL2, respectively. The estimated positions and effects of the detected QTLs affecting HD and PH (denoted IQTLs) are given in Table 5 and Additional files 5 and 6 Tables S5 and S6. The LOD score profiles across the genome in the three populations are shown in Additional file 10: Figure S3. In 8PL, five QTLs were detected for HD, and three were detected for PH (

Comparison of the Three Individual Maps with their Integrated Map
To the best of our knowledge, this study provided the first linkage map from an eight-parent RIL population in rice and the first integrated map constructed from a number of multi-parent populations in rice. Three linkage maps were first constructed from three multi-parent populations, and an integrated map was then built from the three individual maps. Among the three individual maps, the 8PL map had a larger number of bins, shorter genome length, and shorter average bin size than the 4PL1 and 4PL2 maps (Table 1, Figs. 2, 3 and 4). At one locus, the eight-parent population contained twice as many alleles as the four-parent population. The greater number of recombination events and higher genetic diversity in 8PL resulted in more unique locations and higher saturation of markers on linkage maps. The much larger number of markers contained in the integrated map than in the individual maps led to the longer genome of the integrated map. In addition, the integrated map had more bins and a longer length than the three individual maps, with a larger average bin size than that in 8PL but a smaller average bin size than those in 4PL1 and 4PL2. The marker order in the integrated map was not completely consistent with that in the individual maps.
For each linkage map, over ten thousand markers were selected after the first step of filtering, but only half of them were retained after binning. Bin number is much smaller than marker number in linkage mapping populations with limited crossovers, especially when high density markers are used for genotyping. For example, in an actual cowpea population consisting of 305 RILs from an eight-way cross reported by Huynh et al. (2018), 32, 114 SNPs were distributed on 11 chromosomes, but the number of bins was only 1578. One way to improve bin number is to increase population size, by which more recombination events can be observed. Due to experimental and economic limitations, this option may not be applicable in all cases. Integrated maps provide an effective alternative with which to improve bin density, as markers that cannot be distinguished on some maps may have unique positions on other maps (Allen et al. 2017). However, our study failed to show that this is always the case. The integrated map reported in this study was built from three maps of 8PL, 4PL1 and 4PL2, which contained only the markers retained after binning (Table 2). We also tried to construct the integrated map using all markers before binning, i.e., 13,254, 10,825 and 13,479 SNPs in 8PL, 4PL1 and 4PL2. The results showed that the inclusion of all useful markers before binning did not always improve the bin density of the integrated map. For example, for chromosome 12 in the integrated map, the bin number was 674 when using markers before binning, which was smaller than that (677) when using markers after binning. The average bin distance was 0.29 cM, which was larger than that (0.25) observed when using markers after binning.

Comparison of HD and PH QTLs Detected by ICIM and GWAS from the Three Populations
ICIM detected 3, 7 and 9 QTLs for HD and 6, 9 and 10 QTLs for PH in 8PL, 4PL1 and 4PL2 using the three individual maps, respectively. It was not easy to compare the positions of the detected QTLs between maps because of inconsistent marker orders and positions. However, some QTLs had similar physical positions. For example, assuming the length of the confidence interval was 0.5 Mb, the gene Hd3a was detected in both 8PL and 4PL1; Ghd8 was detected in both 8PL and 4PL2; and sd1 was detected in both 8PL and 4PL1 populations.
In a previous study, Shen et al. (2017) applied GWAS in 8PL and detected two QTLs affecting HD and three affecting PH. ICIM detected all these QTLs detected by the GWAS, i.e., qHD6.2, qHD8, qPH1.1, qPH1.2 and qPH12. In addition, ICIM detected one other QTL affecting HD on chromosome 6 and three other QTLs affecting PH on chromosomes 2, 8 and 11, whose PVE were equal to 5.15, 2.13, 1.42 and 1.28%, respectively. In 4PL1, ICIM also detected all QTLs that were detected   by the GWAS, i.e., qHD6 and qPH1. In addition, six other QTLs affecting HD and eight other QTLs affecting PH were identified on various chromosomes. Largest and smallest PVE of these QTL were 11.08 and 1.69%, and average PVE was 3.62%. In 4PL2, two HD and two PH QTLs were identified by the GWAS, one of which was coincident with those identified by ICIM, i.e., qHD8. ICIM detected eight other HD QTLs and ten other PH QTLs. Largest and smallest PVE of these QTL were 34.85 and 0.41%, and average PVE was 4.06%. In summary, GWAS missed some strong signals as well as some weak signals in QTL mapping. The reason may be that GWAS is lack of background control leading to a larger sampling error. The selection of GWAS methods also affects the accuracy of QTL detection. ICIM detected more QTLs, and the detected QTLs explained more phenotypic variance than those detected by GWAS, indicating ICIM has more power in QTL detection than GWAS analysis. For the six cloned genes for HD (i.e. Hd3a, Ghd8 and HD6) and PH (Psd1, sd1 and d61), both GWAS and ICIM detected them for nine times (Additional file 7: Table S7). Hd3a and sd1 was detected by GWAS in 4PL2, but not found by ICIM; HD6 and d61 were detected by ICIM in 4PL2, but not found by GWAS. For the seven times of common genes detected by both methods, the distance of significant SNPs to the gene was smaller by GWAS for five times, and smaller by ICIM for two times. But for most times, the difference between the two methods was small. The reason may be that the most significant SNPs to the gene were deleted by quality control before linkage map construction and QTL mapping. For example, the most significant SNP marker associated with Psd1 in GWAS was not existed in linkage mapping for 8PL.

Comparison of QTLs Detected by Meta-Analysis and QTL Mapping from Integrated Map
Meta-analysis was used to project QTLs detected in the three populations onto the integrated map. In total, 19 MQTLs affecting HD and 25 MQTLs affecting PH were detected by meta-analysis. Some QTLs with similar physical positions across populations were still projected to similar linkage positions on the integrated map. For example, both qPH1.1 in 8PL and qPH1 in 4PL1 were close to the gene sd1 on the rice physical map, corresponding to MqPH1.7 and MqPH1.6 on the integrated map, respectively. The distance between MqPH1.7 and MqPH1.6 was 1.95 cM. Similar results were observed for qHD8 in 8PL and qHD8 in 4PL2. However, some exceptions were detected. For example, both qHD6.2 in 8PL and qHD6 in 4PL1 were close to the gene Hd3a on the physical map, corresponding to MqHD6.2 and MqHD6.4 on the integrated map, respectively. The distance between MqHD6.2 and MqHD6.4 was 31.24 cM. This phenomenon may have been caused by the difference in marker order between linkage and physical maps.
QTL mapping of the three populations based on an integrated map revealed a total of 24 IQTLs affecting HD and 25 IQTLs affecting PH. Compared with the results from the meta-analysis, the numbers of QTLs were similar. Assuming that the length of the confidence interval was 10 cM, 11 QTLs were coincident between MQTLs and IQTLs. For example, MqHD6.1 was located at 10.17 Percentage of phenotypic variance explained cM on chromosome 6, and IqHD6.1 in 8PL was located at 8.10 cM on the same chromosome. The distance between the two QTLs was 2.07 cM. Whether MQTLs or IQTLs are more reliable requires further investigation. Currently, we are developing a combined analysis in which QTL mapping will be conducted by all individuals from each population and a common genetic map constructed using the combined population.

Comparison of Genetic Analysis Using Populations Derived from Different Number of Parents
In this study, the number of QTLs detected in 8PL was lower than that in 4PL1 and 4PL2, along with lower PVE. The reason may be as follows. Degree of freedom of the LOD test statistic in 8PL is higher than that in 4PL. LOD threshold from permutation tests in 8PL is always higher than that in 4PL. Some QTLs detected in 4PL may not reach the threshold level in 8PL, and therefore were not reported when 8PL was used. This may explain why more QTL were detected in 4PL in this study. However, when four parents carry same alleles, this locus cannot be detected in the pureline population derived from the four parents. When other four parents carry different alleles, this locus could be detected in the 8PL population. In this situation, more QTL may be detected in 8PL.Therefore, it is more likely a special case that fewer QTLs were detected in 8PL in this study. Genetic variance at one locus depends on both allele frequencies and genotypic values. PVE is the proportion of genetic variance to phenotypic variance. PVE in 8PL could be higher or lower than that in 4PL in terms of the same locus. For example, given that there are only two genotypic values at one locus, i.e., 10 and 8, and each genotype have equal frequency. In 4PL, if three genotypic values are 10 and one is 8, the genetic variance is equal to 0.75. For the same locus in 8PL, if four genotypic values are 10 and the other four are 8, the genetic variance is equal to 1, which is higher than that in 4PL. If seven genotypic values are 10 and the other one is 8, the genetic variance is equal to 0.4375, which is lower than that in 4PL. In addition, background genetic variance and random error variance could also be different between 4PL and 8PL. Therefore, it is more likely a special case that lower PVE was detected in 8PL in this study. We understand that further studies are needed to make perspective comparison on QTL detection using populations derived from different number of parents.
From the mapping results, it seems that each parent has unique allele and effect in multi-parent populations. In fact, some parents may share same alleles, i.e. the number of unique effect value may be less than the number of parents. However, due to sampling error, the estimated genotypic values are seldom the same even for same alleles; the estimated genotypic values could be similar for genetically different alleles. Therefore, we doubt the exact number of alleles in 4PL or 8PL can be well determined statistically. For convenience, our QTL mapping method assumes different genotypic values for different parents. Despite this, similar effects among parents were observed in many QTLs, which may represent same alleles, for example in Table 3, a 2 and a 3 of qHD6.1, a 1 and a 2 of qHD6.2, a 4 and a 5 of qHD8 and so on. On the other hand, even if we acquired similar estimated effects, it is still hard to determine whether they are the same alleles or not. Instead, biological information is needed to further determine the number of alleles at each QTL, such as DNA sequence and functional analysis.
Compared with bi-parent populations, MAGIC populations contain more alleles and greater genetic diversity. Genetic loci that have no variations in bi-parent populations may have variation in multiple-parent populations because of more alleles are involved. More parents will increase the number of recombination events and improve mapping precision (Huang et al. 2012), which can facilitate QTL detection. However, multiple-parent population is more difficult to develop, because of huge time and labor consumption. More generation of selfing is needed to generate the final pure lines as well. In addition, larger numbers of alleles and marker types at each locus in MAGIC populations complicate the recombination frequency estimation and linkage map construction. Further studies are needed to compare QTL detection using populations from different number of parents.

Conclusions
In summary, individual linkage maps were constructed from three multi-parent rice populations, and an integrated map was built based on the three individual maps. This study produced the first linkage map of an eight-parent RIL population and the first integrated map constructed from multi-parent populations in rice. These maps are necessary for efficient QTL linkage mapping of quantitatively inherited traits. The integrated map supported the comparison and integration of QTLs detected in individual population. QTLs for the traits HD and PH obtained using ICIM based on the three linkage maps were largely coincident with the QTLs detected by a previous GWAS. In addition, ICIM detected more QTLs and explained more phenotypic variance than did the GWAS. The HD and PH QTLs detected by meta-analysis and ICIM based on the integrated map partly overlapped.
Additional file 1: Table S1. Phenotypic mean of the eight parents Additional file 2: Table S2. Information on the individual mapping population data used for linkage maps and QTL analysis