Construction of Consensus Genetic Map With Applications in Gene Mapping of Wheat (Triticum aestivum L.) Using 90K SNP Array

Wheat is one of the most important cereal crops worldwide. A consensus map combines genetic information from multiple populations, providing an effective alternative to improve the genome coverage and marker density. In this study, we constructed a consensus map from three populations of recombinant inbred lines (RILs) of wheat using a 90K single nucleotide polymorphism (SNP) array. Phenotypic data on plant height (PH), spike length (SL), and thousand-kernel weight (TKW) was collected in six, four, and four environments in the three populations, and then used for quantitative trait locus (QTL) mapping. The mapping results obtained using the constructed consensus map were compared with previous results obtained using individual maps and previous studies on other populations. A simulation experiment was also conducted to assess the performance of QTL mapping with the consensus map. The constructed consensus map from the three populations spanned 4558.55 cM in length, with 25,667 SNPs, having high collinearity with physical map and individual maps. Based on the consensus map, 21, 27, and 19 stable QTLs were identified for PH, SL, and TKW, much more than those detected with individual maps. Four PH QTLs and six SL QTLs were likely to be novel. A putative gene called TraesCS4D02G076400 encoding gibberellin-regulated protein was identified to be the candidate gene for one major PH QTL located on 4DS, which may enrich genetic resources in wheat semi-dwarfing breeding. The simulation results indicated that the length of the confidence interval and standard errors of the QTLs detected using the consensus map were much smaller than those detected using individual maps. The consensus map constructed in this study provides the underlying genetic information for systematic mapping, comparison, and clustering of QTL, and gene discovery in wheat genetic study. The QTLs detected in this study had stable effects across environments and can be used to improve the wide adaptation of wheat cultivars through marker-assisted breeding.

Wheat is one of the most important cereal crops worldwide. A consensus map combines genetic information from multiple populations, providing an effective alternative to improve the genome coverage and marker density. In this study, we constructed a consensus map from three populations of recombinant inbred lines (RILs) of wheat using a 90K single nucleotide polymorphism (SNP) array. Phenotypic data on plant height (PH), spike length (SL), and thousand-kernel weight (TKW) was collected in six, four, and four environments in the three populations, and then used for quantitative trait locus (QTL) mapping. The mapping results obtained using the constructed consensus map were compared with previous results obtained using individual maps and previous studies on other populations. A simulation experiment was also conducted to assess the performance of QTL mapping with the consensus map. The constructed consensus map from the three populations spanned 4558.55 cM in length, with 25,667 SNPs, having high collinearity with physical map and individual maps. Based on the consensus map, 21, 27, and 19 stable QTLs were identified for PH, SL, and TKW, much more than those detected with individual maps. Four PH QTLs and six SL QTLs were likely to be novel. A putative gene called TraesCS4D02G076400 encoding gibberellin-regulated protein was identified to be the candidate gene for one major PH QTL located on 4DS, which may enrich genetic resources in wheat semi-dwarfing breeding. The simulation results indicated that the length of the confidence interval and standard errors of the QTLs detected using the consensus map were much smaller than those detected using individual maps. The consensus map constructed in this study provides the underlying genetic information for systematic mapping, comparison, and clustering of QTL, and gene discovery in wheat genetic study. The QTLs detected in this study had stable effects across environments and can be used to improve the wide adaptation of wheat cultivars through marker-assisted breeding.
Keywords: wheat (Triticum aestivum L.), consensus genetic map, QTL mapping, plant height, spike length, thousand-kernel weight INTRODUCTION Wheat (Triticum aestivum L.) is one of the most important cereal crops worldwide, providing about one-fifth of the total calories consumed by humans. Due to limited farmland and the rapid increase in human population, there is an urgent need to accelerate the genetic gain on grain yield through advanced genetic research and breeding activities in wheat. Genetic linkage map construction and quantitative trait locus (QTL) mapping are important areas in genetic research, as they provide fundamental information for gene cloning, marker-assisted breeding, and genome structure studies (Meng et al., 2015;Rasheed et al., 2016).
Linkage mapping approach based on individual populations has become routine in wheat genetic studies to dissect the genetic architecture of complex traits. However, a large number of co-localized markers and low marker density due to a limited genetic variation and a limited number of crossingover events are commonly seen with linkage maps constructed in individual populations. Detected QTLs are usually specific to designated crosses with wide confidence intervals, hindering further genetic research on gene fine-mapping and cloning. Furthermore, linkage mapping in single populations can only identify QTLs with phenotypic variations from specific crosses, and each mapping population can only represent a small number of crossing-over events (Liu and Zeng, 2000). The narrow genetic basis associated with individual crosses and populations reduces both phenotypic and genotypic diversity. One way to solve these problems is to construct a consensus map as the connection across multiple populations.
A consensus genetic map combines genetic information from multiple populations, and therefore provides an effective alternative to improve genome coverage and marker density (Maccaferri et al., 2015;Allen et al., 2017). A higher marker density of the consensus map offers the chance to map more QTLs to narrower intervals and to identify more closely linked markers for the discovery of causal genes and marker-assisted selection (MAS) in breeding. Consensus maps can also be used to validate marker order, characterize genomic diversity, increase the power of genome-wide association studies, and conduct QTL meta-analysis (Cavanagh et al., 2013;Wang et al., 2014;Wingen et al., 2017;Liu et al., 2020).
Some computer tools that can be used for consensus map construction have been developed in the last 20 years, such as BioMercator (Arcade et al., 2004), JoinMap (Van Ooijen, 2006), MergeMap (Wu et al., 2010), MultiPoint (Ronin et al., 2012), and LPmerge (Endelman and Plomion, 2014). Using these tools, consensus maps have been developed for wheat. Somers et al. (2004) reported the first consensus map for wheat based on SSR markers from three doubled haploid (DH) and a recombinant inbred line (RIL) populations. Cavanagh et al. (2013) generated a high-density consensus map from seven populations, consisting of 7,504 single nucleotide polymorphism (SNP) markers. Wang et al. (2014) integrated six bi-parental DH populations to generate a consensus map using 40,267 markers. Liu et al. (2020) developed a consensus map with a total length of 4,080.5 cM containing 47,309 markers based on 21 individual linkage maps and three previously reported consensus maps.
In this study, a consensus genetic map was constructed using three bi-parental populations of RILs in wheat. QTL mapping was then conducted for plant height (PH), spike length (SL), and thousand-kernel weight (TKW) using the constructed consensus map. The mapping results were compared among populations, and with the results obtained using individual maps with the purpose of identifying stable and common QTLs. In addition, a simulation experiment was conducted to demonstrate the advantages of using a consensus map in QTL mapping.

Plant Materials and Phenotyping Experimental Design
The three recombinant inbred line populations used in this study were derived from crosses Doumai × Shi 4185 (denoted as DS, 275 F 2 : 6 RILs), Gaocheng 8901 × Zhoumai 16 (denoted as GZ, 176 F 2 : 6 RILs), and Zhou 8425B × Chinese Spring (denoted as ZC, 245 F 2 : 8 RILs), which had been previously reported by Wen et al. (2017). Population DS and its parental lines were planted at Shunyi (Beijing, China) and Shijiazhuang (Hebei, China) in 2012-2013, 2013 cropping seasons. Population GZ and its parental lines were planted at Anyang (Henan, China) and Suixi (Anhui, China) in 2012-2013 and 2013-2014 cropping seasons. Population ZC and its parental lines were planted at Zhengzhou and Zhoukou (Henan, China) in 2012-2013 and 2013-2014 cropping seasons. Randomized complete block designs with three replications were used in field trials. Each plot had three rows with 1.5 m in length and 0.2 m apart between rows. About 50 seeds were sown in each row. Field management was performed according to local practices.
Plant height was recorded as the average height based on 10 representative plants, measured from the base of the stem to the top of the spike excluding awns at the late grain-filling stage. SL was recorded as the average length of 20 representative spikes in populations DS and GZ, and five representative spikes in population ZC, measured from the base of the spike to the top of the spike excluding awns. TKW was evaluated by weighing three random samples of 500 kernels from each plot after harvest.

Genotyping and Marker Quality Control
Deoxyribonucleic acid was extracted from leaves of 15-day-old seedlings according to the cetyltrimethyl ammonium bromide (CTAB) protocol (Sharp et al., 1988). The populations were genotyped by the 90K wheat Infinium iSelect SNP array  at CapitalBio Corporation (http://www.capitalbio. com) in Beijing, China. Quality control of the genotypic data has been previously described in Wen et al. (2017), and described here briefly and as follows. First, heterozygous marker types were set as missing values. Then, markers with more than 10% of missing values were deleted. Finally, SNPs with minor allelic frequency lower than 0.3 were filtered out. The three individual linkage maps based on these markers were reported by Wen et al. (2017). SNPs on the three maps were used for consensus map construction. The R package VennDiagram (Chen and Boutros, 2011) was used to demonstrate the SNP numbers common among the three individual maps.

Statistical Analysis for Phenotypic Data
Analysis of variance and calculation of broad-sense heritability (H 2 ) from phenotypic data were performed using the AOV function in software QTL IciMapping V4.2 (Meng et al., 2015). Pearson correlation coefficients among traits were calculated using mean phenotypic values across environments.

Consensus Genetic Map Construction
First, markers from the three recombinant inbred line populations were grouped according to their chromosome information in individual maps reported by Wen et al. (2017). Markers that were present on the same chromosome in the three individual maps were treated as anchors. Then, an algorithm called combined linkage analysis (CLA, developed by the group of the authors) was used for consensus map construction. To assure the quality of the map, a limited number of markers were removed manually if they caused serious inconsistency in the marker order between the genetic and physical maps, or excessive expansion of the constructed genetic map. The R package LinkageMapView (Ouellette et al., 2018) was used to visualize the constructed consensus map.
Furthermore, four steps were involved in the CLA algorithm: step 1: derive the theoretical recombination frequencies of pairwise markers in each mapping population; step 2: estimate the recombination frequency between two linked markers and sampling variance of the estimated recombination frequency in each population. In addition to RIL populations, CLA is applicable to many other kinds of bi-parental populations, as described in Meng et al. (2015). For some kinds of mapping populations such as DH and RIL, the likelihood equation on recombination frequency has an explicit solution, so the maximum likelihood estimate can be calculated directly. For other kinds of mapping populations such as F 2 and F 3 , the maximum likelihood estimate cannot be succinctly given. In this situation, either Newton iteration or the expectationmaximization (EM) algorithm has to be adopted in estimating the recombination frequency .
Step 3: estimate the combined recombination frequency using the estimates and their sampling variances from individual populations; reciprocal of sampling variance of the estimated recombination frequency is used as the weight of the corresponding population. Weight is set as zero for those populations where the pair-wise recombination frequency cannot be estimated.
Step 4: construct the consensus linkage map based on the combined estimates of recombination frequencies between markers; a combination of the nearestneighbor algorithm and a two-opt algorithm in solving the traveling salesman problem (TSP) was used in the marker ordering (Zhang et al., 2020a).

Comparison of Marker Orders in the Consensus Map, Physical Map, and Individual Genetic Maps
Spearman rank correlation was used to measure the collinearity of marker orders between the different maps, which was calculated by the R Software. Marker order in each chromosome in the consensus map was compared with the physical map order of the respective chromosome. To acquire the physical positions of the markers, sequences of SNPs were used to BLAST (Basic Local Alignment Search Tool) against the wheat genome IWGSC RefSeq v2.0 (https://urgi.versailles.inra.fr/download/ iwgsc/IWGSC_RefSeq_Assemblies/v2.0/, International Wheat Genome Sequencing Consortium). The E-value threshold in BLAST was set at 10 −10 . The markers were filtered out if their alignment lengths were lower than 80% of the query sequence length or the identities were lower than 0.85. If a marker was assigned to multiple chromosomes by BLAST, the position on the same chromosome as the consensus map was used in collinearity analysis. Marker order comparison was also conducted between the consensus map and individual maps, as well as among the three individual maps. For each comparison, only the common markers on two maps were used in the calculation of collinearity.

QTL Mapping Based on the Consensus Map
Quantitative trait locus mapping was conducted in the individual populations using the consensus map. The inclusive composite interval mapping (ICIM) implemented in the BIP function in QTL IciMapping V4.2 (Meng et al., 2015) was applied on the mean phenotypic values across blocks in each environment and best linear unbiased estimation (BLUE) values across multiple environments. Scanning step was set at 0.2 cM. Probabilities of adding and removing variables in stepwise regression were set at 0.001 and 0.002, respectively. Threshold logarithm of odd (LOD) score was set at 2.5, same as the QTL mapping studies on individual maps from the three populations (Gao et al., 2015;Li et al., 2018).
Quantitative trait loci and quantitative trait locus clusters were named with chromosomal locations, considering all populations together. QTLs detected in the same population were considered to be common if the distance between QTL positions was <20 cM in different environments. QTLs detected in different populations were considered to be common if the genetic and physical positions were close enough. In other words, distance in the linkage map was <20 cM in terms of QTL positions, and distance in the physical map was <25 Mb in terms of the minimum physical distances between flanking makers. In individual populations, QTLs are considered to be stable if they are identified in at least half of tested environments. Stable QTLs for different traits were classified into the same cluster if the minimum distance between the QTL confidence intervals was <15 cM. The shiny Circos tool (Yu et al., 2018) was used to visualize QTL positions on the consensus map. Stable QTLs detected with the consensus map in this study were compared with those detected by ICIM using individual maps (Gao et al., 2015;Li et al., 2018), according to physical and genetic positions of the flanking markers.

Genetic Models Used in Simulation
A simulation study was conducted to compare the QTL mapping results from the individual and consensus maps. We assumed that a chromosome has a length of 100 cM and contains 101 evenly distributed markers. Considering that the consensus map always has more markers than each individual map, we assume that the consensus map contained all the 101 markers, but that the individual map only contained half of them, i.e., 51 evenly distributed markers with marker density at 2 cM. Three QTL  (Meng et al., 2015). The consensus map with 101 markers and the predefined QTLs were used to generate the simulated populations. Both the consensus and individual maps were used in QTL mapping. For QTL mapping using individual maps, genotypic data of the 51 markers were used. For QTL mapping using the consensus map, genotypic data of the 51 markers were the same as those in individual maps, but the other 50 markers only present in the consensus map were set as missing values. For the ICIM QTL mapping method on simulated populations, the scanning step was set at 0.1 cM and the threshold LOD score was set at 2.5. Probabilities for entering and removing variables in the stepwise regression were set at 0.001 and 0.002, respectively. QTL detection power was estimated according to a support interval of 5 cM centered at the position of true QTL. If multiple peaks occurred within the support interval, only the highest one was counted. QTLs identified out of the support interval were regarded as false positives (Li et al., 2012). The other parameters were set as default values.

General Information on Both Genotypic and Phenotypic Data
There were 10,986 markers on the linkage map constructed in population DS, 11,819 markers in population GZ, and 14,862 markers in population ZC. Populations DS and GZ shared 4,208 common markers; DS and ZC shared 4,420 common markers; GZ and ZC shared 5,183 common markers; the three populations had 1,880 markers in common (Supplementary Figure 1). A total of 25,736 unique markers on the three individual maps were used for consensus map construction. Phenotypic means and heritability of the three traits are shown in Table 1 Figure 2). For SL, Doumai was longer than Shi 4185 in four environments, almost equal in one environment, and shorter in the other one environment; Zhoumai16 was longer than Gaocheng 8901 in three environments, and shorter in the other environment; Chinese Spring was always longer than Zhou 8425B (Supplementary Figure 3). For TKW, Doumai was always higher than Shi 4185; Zhoumai 16 was always higher than Gaocheng 8901; Zhou 8425B was always higher than Chinese Spring (Supplementary Figure 4). The three traits were continuously distributed in the three populations, similar to and typical in most QTL mapping studies. Much wider ranges were observed in the progenies in comparison with their parents, except for TKW in two environments in population ZC (Supplementary Figures 2-4). Heritability in the broad sense, based on the replicated means, was quite high for the three traits, ranging from 0.91 to 0.99 (Table 1), while heritability based on the plot level ranged from 0.59 to 0.91. Correlation coefficients between traits in the three populations are given in Supplementary Table 2. At a significance level of 0.01, PH was positively correlated with TKW in all three populations. SL showed a positive correlation  with both PH and TKW in population DS. Other correlations were non-significant.

Characteristics of the Constructed Consensus Map
Of the 25,736 unique SNPs on the three individual linkage maps, 25,667 were assigned to the consensus map, resulting in 21 linkage groups corresponding to the 21 chromosomes in hexaploid wheat (Figure 1). General information on the consensus map is provided in Table 2, and positions of all the markers on both the genetic and physical maps are given in Supplementary Table 3. The consensus map spanned 4,558.55 cM in length, and the number of unique map positions (denoted as bins) was equal to 3,979. Lengths of the A, B, and D genomes were 1,622.47, 1,581.57, and 1,354.52 cM, respectively ( Table 2). Chromosome 4D was the shortest, with a length of 168.06 cM, and had the least number of markers (i.e., 106) and the least number of bins (i.e., 53). Chromosome 2B was the longest with a length of 290.78 cM, and had the second largest number of markers (i.e. 2,431) and the second largest number of bins (i.e., 324). There were 18 gaps longer than 15 cM on the consensus map, 16 of which were located in the D genome (Supplementary Table 3). Average distance between adjacent bins was equal to 1.15 cM. The single nucleotide polymorphism markers (SNPs) number was similar in the A and B genomes, i.e., 10,285 and 12,755 SNPs, but the number was much lower in the D genome, i.e., 2,627 SNPs ( Table 2). In comparison with the A and B genomes, the D genome was shorter and contained much fewer markers and bins, and more gaps, indicating that fewer crossing-over events happened on the D genome, which was also observed in the three individual maps. Although the marker number and bin number in the D genome were significantly lower than those in the A and B genomes, results from BLAST indicated that the constructed consensus map still had nearly complete coverage for chromosomes in the D genome.
Marker orders on the consensus map and physical map had high collinearity, with an average Spearman rank correlation coefficient of 0.95 across the 21 chromosomes (Table 2, Figure 2). Rank correlation coefficients were higher than 0.94 for all the chromosomes except 2B and 3D. The lower coefficient observed on 3D may be partly due to the much-reduced bin number when many markers were clustered in bins. Collinearity analysis between the consensus and physical maps also revealed that markers in large physical region around the centromeres of chromosomes tended to be clustered in a short genetic interval on consensus genetic map (Figure 2), indicating a much stronger recombination suppression occurred around the centromere than did that the distal regions. Wen et al. (2017) reported three linkage maps from three populations constructed with QTL IciMapping V4.0 (Meng et al., 2015), JoinMap 4.0 (Stam, 1993), and MapDisto 1.7 (Lorieux, 2012). Two of them had 21 linkage groups, and one had 31 linkage groups. The consensus map constructed in this study had 21 linkage groups, corresponding to the 21 chromosomes in hexaploid wheat. The marker and bin numbers on the consensus map were 1.73 and 1.15 times higher than the largest marker and bin numbers on the three individual maps. The length of the consensus map was 1.44 times longer than that of the longest individual map. Longer chromosomes on the individual maps also tended to be longer on the consensus map. For example, the two longest chromosomes on the consensus map, i.e., 2B and 5B, ranked first and third in mapping length in each of the three individual maps.

Comparison of the Consensus Map With the Three Individual Maps
There were 616 markers with inconsistent chromosomes on the individual maps, but the inconsistent chromosomes for each marker were finalized to one unique chromosome on the consensus map (Supplementary Table 4). Among these markers, 540 were mapped to single chromosomes that they were located on the individual maps. For example, marker wsnp_Ex_c200_391015 was located on chromosomes 7A and 1A on individual maps of populations GZ and ZC, respectively, which was finalized on chromosome 1A on the consensus map. Forty-nine markers were mapped to one of the homeologous chromosomes. For example, marker Tdurum_contig28665_150 was located on chromosomes 1D, 1D, and 2A in populations DS, GZ, and ZC, respectively, and was finalized on chromosome 1A, a homeologous chromosome of 1D. Twenty-seven markers were mapped to neither the same chromosome nor homeologous chromosomes. For example, marker tplb0024a09_2369 was located on chromosomes 7D and 4A in populations DS and ZC, respectively, and was finalized on chromosome 2B on the consensus map (Supplementary Table 4).
The markers showed high collinearity across chromosomes between the consensus and individual maps, and the average Spearman rank correlation coefficient was similar to those between the individual maps (Supplementary Table 5). Fewer inconsistencies in orders between the consensus and individual maps were observed for closely linked markers.

QTLs for TKW Detected From the Consensus Map and Comparison With Those From the Individual Maps
Using the consensus map, a total of 53 QTLs were detected for TKW (Supplementary Table 6), among which nine, three, and eight were stable in populations DS, GZ, and ZC, respectively (Figure 3, Table 5). qTKW-4B-2 was repeatedly identified in populations DS and GZ with LOD scores ranging from 3.08 to 49.22, explaining 7.57-36.51% of the phenotypic variance. qTKW-4B-2 had the largest LOD score, PVE and additive effect across environments in population DS. This QTL was co-localized with qPH-4B-1, corresponding to the dwarfing gene Rht-B1. All stable QTLs detected with the individual maps were also stable when detected with the consensus map (Supplementary Table 7, Supplementary Figure 7). There were other three stable TKW QTLs identified using the consensus map, i.e., qTKW-1B-2, qTKW-2D-2, and qTKW-6B-3. qTKW-1B-2 was mapped on chromosome 1B at the interval of 588.

QTL Clusters for the Three Traits
As far as the stable QTLs across environments were concerned, 11 QTL clusters were identified and distributed on nine chromosomes (Supplementary Table 8), six of which affected two traits (i.e., qClu-2D, qClu-4A-1, qClu-4A-2, qClu-4D, qClu-6A, and qClu-6B), and five affected all the three traits (i.e., qClu-3A-1, qClu-4B, qClu-5A-1, qClu-5A-2, and qClu-7A). Eight clusters affected traits PH and SL. Among them, three clusters contained both PH and SL QTLs in population DS; one cluster contained both PH and SL QTLs in population ZC, and one cluster contained the closely linked PH and SL QTLs in populations DS and GZ. Each of the five clusters either increased or decreased both traits simultaneously. Genomic regions containing the stable QTLs for the three traits were located on chromosomes 3A, 4B, 5A, and 7A. The cluster on 4B was close to the Green Revolution gene Rht-B1. In cluster qClu-5A-1, QTLs affecting the three traits were consistently identified in populations DS and GZ, either increasing or decreasing the three traits simultaneously.

Potential Applications of the Detected QTLs in Wheat Breeding
To explore the potential applications of the detected QTLs in wheat breeding, QTL genotypes and genotypic values of each RIL in the three populations were predicted on the three traits with stable QTLs identified using BLUE values across environments (Supplementary Tables 9-11). For convenience, for the two alleles at each QTL, one is called positive and the other one is called negative. Parental sources of the two alleles can be determined from the sign of the estimated additive effect of the QTL. Due to the varied objectives on different traits in breeding, it should be noted that the positive allele is not always favored and that the negative allele is not always un-favored. For PH, nine, eight, and eight stable QTLs were used for prediction in populations DS, GZ, and ZC, respectively. The 10 highest RILs possessed at least eight, seven, and seven positive alleles in the three populations, respectively, whereas the 10 lowest RILs had no more than two positive alleles (Supplementary Table 9). For SL, 14, 7, and 4 stable QTLs were used for prediction. The 10 highest RILs possessed at least nine, seven, and four positive alleles in the three populations, whereas the 10 lowest RILs had no more than four positive alleles in population DS, no positive allele in population GZ, and no more than 1 positive allele in population ZC (Supplementary Table 10). For TKW, seven, three, and six stable QTLs were used for prediction. The 10 highest RILs possessed at least six, three, and five positive alleles in the three populations, whereas the 10 lowest RILs had no more than 1 positive allele (Supplementary Table 11 Recombinant inbred lines with the predicted genotypic values on PH, SL, and TKW can serve for the choice of target genotypes meeting different breeding objectives, such as wheat cultivars with medium plant height, large spike length, and medium to high kernel weight. Given one target genotype, the predicted allelic combination of RILs can serve for the prediction of cross performance and the selection of suitable parental lines through simulation or other genomic prediction approaches .

QTL Mapping in Simulated Populations
In 1,000 simulated populations, the estimated QTL positions and effects using the individual and consensus maps are shown in Table 6. With the increase in heritability, QTL detection powers were increased and the false discovery rate (FDR) was decreased in the three models using either the individual or consensus maps. Approximately unbiased estimation of QTL positions and effects was obtained for each defined model and heritability level. The confidence intervals of QTLs detected from the consensus map were much narrower, and the associated standard errors were much smaller than those from individual maps. Detection power was much lower for QTLs in linkage models II and III than that in the unlinked model I at the same heritability levels for both the individual and consensus maps. FDR was much higher in models II and III than in model I, indicating the complexity and difficulty in dissecting linked QTLs in genetic studies.

Computer Tools in Consensus Map Construction
Two strategies have been adopted for consensus map construction in previous studies (Endelman and Plomion, 2014). The first one is based on the raw data of multiple mapping populations, and has been implemented in software MultiPoint (Ronin et al., 2012) and JoinMap (Van Ooijen, 2006). The second one is based on individual linkage maps previously constructed, and has been implemented in software BioMercator (Arcade et al., 2004), MergeMap (Wu et al., 2010), LPmerge (Endelman and Plomion, 2014), and QTL IciMapping (Meng et al., 2015). The first strategy is usually time-consuming when dealing with a large number of markers (Wu et al., 2010), which has drastically restricted the use of a large number of markers in the consensus map. The second strategy highly depends on the quality of individual maps and sometimes may result in maps with unreasonable length (Cavanagh et al., 2013;Wang et al., 2014;Wingen et al., 2017).
With the development of high-throughput sequencing technology, markers that can be used in genotyping mapping populations are growing rapidly. A large amount of markers brings a great challenge to consensus map construction, especially when raw genotypic data are used. The two raw databased software packages mentioned above cannot deal with such a large number of markers used in this study. For example, both packages cannot generate a consensus map for chromosome 5B, which harbored 929, 1,406, and 1,508 SNPs in populations DS, GZ, and ZC, respectively. Map-based method only utilizes marker distances between adjacent markers, which may result in an inaccurate estimation of recombination frequency between markers especially when the order of markers changes on the consensus map. The CLA algorithm is a raw data-based method used in this study to deal with a large amount of markers. The combined recombination frequency between any pair of markers was calculated from the estimates in individual mapping populations. The estimated recombination frequencies are recorded in computer memory. Therefore, time can be greatly saved in computing.

Quality of the Consensus Map
The great number of markers and bins contained in the consensus map provided higher saturation of markers and better genome coverage, and expanded the length of the map. Previous studies have shown that increased recombination events and map resolution with an increased number of markers and density could contribute to longer map length (Ferreira et al., 2006;Wingen et al., 2017). The longer map length may also suffer from chromosomal structure differences in different mapping populations and the ordering algorithm used. Compared with the A and B genomes, the D genome had fewer unique markers, larger gaps, and shorter map length, which have been previously reported in both consensus and individual maps in wheat Li et al., 2015;Guan et al., 2018). Collinearity was high between the genetic and physical positions. Marker order on the consensus and physical maps was highly correlated at the genome-wide level, but lower collinearity was sometimes observed in some chromosomal regions, which was also reported previously (Wingen et al., 2017). Of the 19,320 SNPs on the consensus map that had physical positions, on average there were 55.17% SNPs arranged in the same order as those on the corresponding chromosomes of the physical maps, ranging from 36.3 on chromosome 6A to 75.68% on chromosome 4D ( Table 2). A higher proportion of the completely consistent marker order was found in the D genome (63.95%) than those in the A genome (52.32%) and the B genome (49.22%), which may be explained by the lower recombination on the D genome. The lower recombination events on the D genome contributed to lower sequence variability and had a weaker influence on the decay of syntenic block size. Some chromosomal structural variations were observed on the consensus map, such as intra-chromosomal translocation and inversion. For example, inversion happened around 22-25 Mb on chromosome 1A, and translocation occurred between regions around 88-93 and 106-109 Mb on chromosome 2A. The collinearity between marker orders in genetic and physical maps is often disturbed by the macrostructural variations in wheat, especially for consensus maps that are constructed from multiple populations. Local disorder of markers could also be caused by the variation of gene order in parents and genotyping errors.
The distribution of meiotic recombination events showed that recombination happened much more frequently in the distal chromosomal regions, and that recombination tended to be suppressed near the centromeres, which was consistent with previous studies [Sourdille et al., 2004;International Wheat Genome Sequencing Consortium (IWGSC), 2018]. Collinearity analysis also showed that some markers might have conservative orders across populations, since their relative orders were consistent on the physical and genetic maps. Comparative analysis among the consensus, physical, and individual maps indicated the reliability of the consensus map constructed with the CLA algorithm.

Comparison of the Detected QTLs With Studies on Other Mapping Populations
In this study, eight stable PH QTLs were detected with the consensus map but not with the individual maps ( Table 7). Guan et al. (2018) reported a PH QTL on chromosome 4D at the physical interval of 37. 05-62.94 Mb, and Ren et al. (2021) reported a PH QTL on the same chromosome at the physical interval of 47.44-67.64Mb. qPH-4D-2 (chr4D:32.97-65.01 Mb) was overlapped with the loci reported by Guan et al. (2018) and Ren et al. (2021). qPH-6A-1 was located within the physical region as reported by Zanke et al. (2014). qPH-6A-2 was mapped on chromosome 6A at the interval of 610. 97-613.55 Mb. Similarly, Pang et al. (2020) detected a PH QTL on chromosome 6A at the interval of 609.3-609.9 Mb (IWGSC RefSeq v1.0). qPH-6D-2 was located at the same marker interval of a PH QTL that was first reported and validated to be stable in two wheat populations by Wang et al. (2020). To the best knowledge of the authors, stable QTLs qPH-2D-1, qPH-2D-3, qPH-3B-2, and qPH-7A identified in this study were likely to be novel for PH. The increased marker density in the consensus map contributed to the detection of these novel QTLs.
For the three traits, a total of 21 QTLs were identified using the consensus map but not the individual maps. Among them, 11 QTLs are consistent with those from previous studies on other mapping populations, and 10 QTLs are likely to be novel. Most of the 11 QTLs were first reported in recent years using high-density linkage maps, indicating that the increase in marker density improved the power of QTL detection. For the novel QTLs, six of them that control PH or SL were included in the cluster that harbored closely linked PH and SL QTLs (Supplementary Table 8). The PH of the wheat plant is equal to SL plus the lengths of all internodes above the ground. Theoretically, loci associated with SL may affect PH as well, which has been validated by some studies (Buerstmayr et al., 2011;Lv et al., 2014;Xu et al., 2014;Jahani et al., 2019;Chen et al., 2020). Furthermore, four novel SL QTLs were close to PH QTLs that have been reported using individual maps or other independent studies, indicating the reliability of the novel QTLs on SL or PH. Gene TaERF8 was identified to be associated with PH and yield in wheat, and has been cloned from the wheat cultivar Chinese Spring (Zhang et al., 2020b), one parental line of population ZC. TaERF8-2D (chr2D: 368.21 Mb) was located in the flanking marker interval of qPH-2D-1, which was stably detected in population ZC in the three tested environments and in population DS in two tested environments. TaERF8-2D may be a candidate gene for qPH-2D-1. Annotations of gene functions were also performed for these novel QTLs based on the wheat reference sequence annotation database (IWGSC Annotation v1.1) as listed in Supplementary Table 12. The annotation information will facilitate the future fine mapping, map-based cloning, and functional analysis of the novel QTLs identified in this study.

Relationship Between QTLs for Phenotypically Correlated Traits PH and SL
Plant height is an important agronomic trait highly related to lodging resistance and harvest index in wheat. SL is highly related to grain yield by affecting kernel number and spike morphology (Donmez et al., 2001). Plants with suitable PH and larger spike are desirable in wheat breeding. Nine of the 21 stable PH QTLs were close to the stable SL QTLs (Supplementary Table 8), contributing to the genetic correlation between the two traits. PH and SL were positively correlated by phenotypic analysis in population DS, but the correlation was non-significant in the other two populations. In this study, closely linked PH and SL QTLs identified in the same population always had genetic effects at the same directions on both traits. Similar instances have been reported in previous studies (Buerstmayr et al., 2011;Lv et al., 2014;Xu et al., 2014;Jahani et al., 2019;Chen et al., 2020). Considering that some QTLs for SL may also affect PH, we speculated that the closely linked PH and SL QTLs are more likely to be the same genetic loci and have the same effect directions.  E1, 2012E1, -2013E2, 2012E2, -2013E3, 2013E3, -2014E4, 2013E4, -2014E5, 2014E5, -2015E6, 2014E6, -2015E7, 2012E7, -2013E8, 2012E8, -2013E9, 2013E9, -2014E10, 2013E10, -2014E11, Zhoukou2013;E12, Zhengzhou 2013;E13, Zhoukou 2014;E14, Zhengzhou 2014;B, best linear unbiased estimation. However, whether the closely linked QTLs on PH and SL belong to the same chromosomal loci with pleiotropic effects or different closely-linked loci needs further investigation and is beyond the scope of this study.

Further Analysis for a Major PH QTL Located on Chromosome 4DS
For plant height, only one QTL was detected on chromosome 4DS using the individual maps in populations GZ and ZC, but two stable QTLs, i.e., qPH-4D-1 and qPH-4D-2, were identified using the consensus map in the same two populations, which were linked in the coupling phase (Table 3,  Supplementary Table 7). The BLAST results indicated that qPH-4D-1 was co-localized with the dwarfing gene Rht-D1. qPH-4D-2 explained 8.39-10.9 and 5.03-12.03% of the phenotypic variance across environments in populations GZ and ZC, respectively. The alleles decreasing PH were from parents Zhoumai16 in population GZ and Zhou 8425B in population ZC. Guan et al. (2018) reported two QTLs that were also linked in the coupling phase and located in similar positions as qPH-4D-1 and qPH-4D-2. qPH-4D-2 was detected in four environments and with BLUE values across eight environments in Guan et al. (2018). In addition, qPH-4D-2 was closely linked with marker wsnp_Ex_c683_1341113, which was also observed in Guan et al. (2018). As reported by Ren et al. (2021), qPH-4D-2 was detected in the similar position between SNPs AX-89692818 and AX-109606880 across environments. Therefore, it is highly possible that qPH-4D-2 is a novel semidwarfing gene. The common marker wsnp_Ex_c683_1341113 was located at about 54.4 Mb on chromosome 4D (IWGSC RefSeq v1.0; IWGSC, 2018). A high confidence putative gene, TraesCS4D02G076400 (50,888,889,461 bp), is located around the marker and in the confidence interval of qPH-4D-2, with the annotation of encoding gibberellin regulated protein (IWGSC RefSeq v1.1 annotation; IWGSC, 2018). Gibberellin is an essential endogenous regulator in plant growth. The well-known dwarfing genes Rht-B1b and Rht-D1b regulate DELLA proteins in gibberellin signaling to reduce the response to gibberellin (Peng et al., 1999). The gibberellin-sensitive gene Rht8 was also widely used in regulating PH in wheat (Gasperini et al., 2012). Gene TraesCS4D02G076400 in wheat was annotated to gene GAST1 (UniProtKB/TrEMBL; Acc:C8C4P9), first reported in tomato to encode the gibberellins-stimulated transcript (Shi et al., 1992). GAST1 belongs to the gibberellic acidstimulated Arabidopsis (GASA) family, which plays important roles in plant growth and development, such as stem growth, plant height, and grain length, width, and weight (de la Fuente et al., 2006;Nahirñak et al., 2012a,b;Shi et al., 2020). Furthermore, qPH-4D-2 was detected in two populations in this study, one from the cross between Zhou 8425B and Chinese Spring. TraesCS4D02G076400 had high RNA expression levels in Chinese Spring in different tissues and development stages (expVIP, http://www.wheat-expression.com/). Therefore, TraesCS4D02G076400 is likely to be the candidate gene for qPH-4D-2. PH is a crucial trait for morphogenesis and grain yield in wheat. The newly discovered PH QTL on chromosome 4DS in this study may enrich the genetic resources in breeding for semidwarfing wheat. Reasons that qPH-4D-2 was not identified by the individual could be the short distance between the QTL and Rht-D1, and the lower marker density around the two QTLs in individual mapping populations.

Advantages of Using Consensus Map in QTL Mapping
Due to the limited number of crossing-overs and limited genetic variation in individual populations, linkage maps constructed from individual mapping populations usually have a large number of co-localized markers and low marker density. A consensus map combines the genetic information included in multiple populations and provides a better genomic coverage with higher marker density (Maccaferri et al., 2015;Allen et al., 2017). A consensus map of higher density offers the chance to map QTLs to narrower chromosomal intervals, which will facilitate the discovery of causal genes and the identification of closely linked markers for MAS. Simulation results conducted in this study confirmed that the use of a consensus map with higher marker density reduced the confidence interval of detected QTLs.
Even for the same trait, QTLs detected in different populations using their own genetic maps sometimes are hardly compared and synthesized, because of the unshared markers and variations in the genetic background (Sukumaran et al., 2015). Comparisons on QTL positions estimated from different populations are usually conducted by anchoring the linked markers to the genome assembly. However, genome sequences usually have wide variations between parental varieties, and the anchor information to the genome sequence may not be completely accurate. A consensus map provides the direct comparison for QTLs detected from different populations, which is important, particularly in species lacking a completely sequenced reference genome. In this study, we demonstrated that QTL mapping using a consensus map can better identify common and stable QTLs across populations and environments. For example, Rht-B1 and Rht-D1, which had been cloned, were the two genes reducing plant height in wheat (Peng et al., 1999). Each of them was located almost in the same position in two populations on the consensus map. qPH-5A-2, qSL-2D-2, qSL-5A-2, and qTKW-4B-2 were detected in populations DS and GZ; qPH-2B-2 was detected in populations DS and ZC; qPH-4B-1, qPH-4D-1, qPH-4D-2, and qSL-2D-1 were identified in populations GZ and ZC; qSL-6B-4 was detected in all the three populations. The common QTLs identified in multiple populations reflected the stable genetic effects of QTLs in different genetic backgrounds, which might be more valuable in breeding.
The genetic relationship among PH and SL QTLs as observed in this study, showed that QTL mapping using the consensus map can also facilitate the comparison across the correlated traits, and therefore provide the opportunity to understand the genetic correlation between phenotypically correlated traits and identify the QTL-rich genomic regions. Moreover, the consensus map also provides the chance to detect common QTLs with smaller effects occurring in different populations.
Further studies may still be needed to determine the key factors affecting the accuracy of consensus map construction and subsequent QTL mapping, such as proportion of common markers shared by multiple mapping populations, inconsistency degree of marker orders in individual populations, populationspecific recombination frequencies, and the optimum algorithm used to construct the consensus map. In addition to bi-parental populations, as have been used in this study, multi-parental populations have been developed in recent years in crops together with suitable genetic analysis methods (Gardner et al., 2016;Zhang et al., 2017Zhang et al., , 2019Shi et al., 2019;Qu et al., 2020). In theory, a consensus map can also be constructed by combining a number of bi-parental and multi-parental populations, when common markers are shared by these populations.
In conclusion, the consensus map constructed for this study allows for systematic QTL mapping studies, and comparison and clustering of mapping results in wheat genetic studies. The QTL mapping based on the consensus map resulted in higher accuracy, narrower confidence interval, and a larger QTL number. The stable QTLs across tested environments and mapping populations, and the predicted QTL genotypes and genotypic values can be used to select wheat cultivars with suitable PH, large SL, and medium to high kernel weight. SNPs closely linked with these stable QTLs can be used to select suitable genetic materials and make suitable crosses in wheat breeding programs. SNPs closely linked to traits can also be converted into Kompetitive allele-specific PCR (KASP) markers (Kaur et al., 2021) and then used for large-scale genotyping to screen desirable individuals in segregating breeding populations.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
LZ and JW conceived and designed the research. PQ and LZ conducted data analysis. FG, WW, JL, and XX developed the populations, performed SNP genotyping, and conducted field trials. PQ, JW, and LZ wrote, drafted, and revised the manuscript. XX and HP provided guidance on data analysis and revised the manuscript. All the authors read and approved the final version of the manuscript for publication.

FUNDING
This study was supported by the National Natural Science Foundation of China (Project No. 31861143003).