Construction of a high density genetic map and QTL mapping of PEG-induced drought tolerance at the early seedling stage in sesame using whole genome re-sequencing

Junchao Liang Crop Research Insitute, Jiangxi Academy of Agricultural Sciences, Nanchang, China https://orcid.org/0000-0003-0635-5243 Yanying Ye Crop Research Institute, Jiangxi Academy of Agricultural Sciences, Nanchang, China Xiaowen Yan Crop Research Institute, Jiangxi Academy of Agricultural Sciences, Nanchang, China Tingxian Yan Crop Research Institute, Jiangxi Academy of Agricultural Sciences, Nanchang, China Yueliang Rao Crop Research Institute, Jiangxi Academy of Agricultural Sciences, Nanchang, China Hongying Zhou Nanchang Institute of Technology Jian Sun (  whsunjian@aliyun.com ) Crop Research Institue, Jiangxi Academy of Agricultural Sciences, Nanchang, China https://orcid.org/0000-0002-7443-179X Meiwang Le Crop Research Institue, Jiangxi Academy of Agricultural Sciences, Nanchang, China


Abstract
Background Improvement in sesame (Sesamum indicum L.) drought tolerance at seedling stage is important for yield stability. Genetic approaches combing with conventional breeding is the most effective way to develop drought-tolerant cultivars. So far, only a few studies have been reported to reveal gene/ quantitative trait loci (QTL) controlling drought tolerance in sesame. To identify the genomic regions associated with drought tolerance, we constructed a high-density genetic map using a recombinant inbred line (RIL) population through whole genome re-sequencing (WGRS) technique. QTLs contributing to three seedling traits were identi ed under both non-stress and water stress conditions.

Results
Three drought tolerance related traits and their relative values (the ratio of value under stress to value under control condition), including seedling weight (SW), shoot length (SL) and root length (RL), were evaluated under control and PEG-induced osmotic conditions at seedling stage in a RIL population derived from cross of Zhushanbai (ZSB) and Jinhuangma (JHM). Signi cant variation and high broad sense heritability were observed for all traits except SW under stress condition in the population. With this population, a highdensity linkage map with 1354 bin markers was constructed through WGRS strategy. Composite interval mapping analysis was performed for all the traits as well as their relative phenotypic data. A total of 34 QTLs were detected for these traits under both conditions and their relative values, and 13 stable QTLs associated with seven traits were revealed in two independent experiments, explaining on average, 4.95-16.26% of phenotypic variation for each QTL. Four of them contributed more than 10% of phenotypic variation. One region on chromosome 12 contained two major QTLs related to RL under osmotic condition and relative RL. Seven candidate genes underlying major QTLs for drought tolerance were identi ed according to gene descriptions and variations between parents.

Conclusion
The current study reports the rst QTL mapping of drought tolerance related traits through a RIL population and rst QTL detection of root related trait (root length) in sesame. These ndings will provide new genetic resources for molecular improvement of drought tolerance and candidate gene identi cation in sesame.

Background
Drought stress usually refers to a water shortage which causes dramatic morphological, biochemical, physiological, and molecular changes [1]. These changes would severely affect plant growth and crop yield stability. In past four decades, it is estimated that the drought has caused a cereal loss of 1820 million Mg [2]. Meanwhile, the continuous global warming and climate change will probably increase the frequency of drought. Sesame (Sesamum indicum L., 2n = 26) is one of the most important oilseed crops in the world. Sesame seeds are known to be rich in protein, vitamins and special antioxidants, such as sesamin and sesamolin, which make sesame become a very healthy food favored by consumers. Comparison with most of other oilseed crops, sesame is considered as a resilient crop that is more tolerant to drought stress. However, sesame growth and productivity are vulnerable to severe drought stress, especially in the arid and semi-arid areas. Sesame belongs to shallow root plants and is very sensitive to drought stress during the germination and owering stages [3]. Improvement of drought tolerance at these two stages of sesame is very important for yield stability.
Drought tolerance is a complex quantitative trait in plants. QTL mapping and genome-wide association study (GWAS) have been widely used for genetic analysis of drought tolerance related traits in many plants, such as rice [4,5], wheat [6][7][8], cotton [9,10] and Brassica napus [11]. These traits including root length, coleoptile length and shoot length for seedling stage, and yield-related traits for owering stage. A solution of polyethylene glycol (PEG) is frequently used to simulate drought stress through treating the seeds with the PEG solution for days, especially in germination or seedling stress experiments [4][5][6]12,13].
In sesame, several studies have reported that many traits including germination rate, seedling growth, shoot length, root length and yield related traits could be affect by drought stress [14][15][16]. Mensah et al. [14] used varying PEG concentrations to simulate drought effect on germination of sesame and found that higher osmotic conditions (0.25-0.50 MPa) signi cantly reduced the germination rate, radical and shoot development, but lower osmotic tensions (0.0625 MPa) could enhance root growth. Only a few studies have been conducted yet on genetic analysis of drought tolerance. Li et al. [17] con rmed 15% PEG 6000 as a suitable concentration for examining drought tolerance in sesame germplasms, and performed a GWAS of stress tolerance indexes related to NaCl-salt and PEG-drought at the germination stage with 490 sesame accessions, and identi ed nine and 15 QTLs for drought and salt stresses, respectively. Ten stable QTLs were also identi ed for ve drought related traits at the sesame owering stage through GWAS [18]. By using gene association study, gene expression and transgenic experiments, a candidate gene SiSAM was identi ed to confer drought tolerance by modulating polyamine levels and ROS homeostasis in sesame [18]. So far, there is no study on QTL mapping for drought tolerance traits using bi-parental population in sesame.
The recent rapid development of high throughput sequencing makes it easier to construct high-density genetic map with numerous single-nucleotide polymorphism (SNP) markers. SNPs are the most abundant form of genetic variation throughout the genomes and are ideal genetic markers for genetic and breeding applications. Recent released sesame reference genome has greatly helped in SNPs identi cation in sesame through multiple next generation sequencing strategies, including genotyping by sequencing (GBS), reduced representation sequencing (RAD) and whole genome re-sequencing. These methods have been successfully utilized for high resolution mapping of QTLs in sesame [19][20][21][22][23].
To further explore the genetic foundation of sesame drought tolerance, in the present study, we developed a high-density genetic map through WGRS of 180 RILs generated by two sesame landraces Jinhuangma (JHM) and Zhushanbai (ZSB) from China. By performing PEG stress experiments, QTLs for fresh seedling weight, shoot length and root length were identi ed under control and stress conditions. These can provide a valuable contribution to understand the genetic basis of drought tolerance and facilitate marker-assisted breeding for stress tolerance in sesame.

Results
Whole genome re-sequencing and genotyping Using whole genome re-sequencing approach, we generated over 1.1 billion reads from two parents and 180 RILs, with an average of 6.3 million reads per RIL, providing an average read depth of 4.18×. For the two parents JHM and ZSB, ~14.6 million and ~13.9 million reads were obtained respectively. All the clean reads were mapped to the sesame reference genome Zhongzhi No. 13 [20] using BWA. After ltering, a total of 466,911 high-quality SNPs and 72,981 InDels (Insertion and deletions) were identi ed among the RILs.  Table S1).

High density genetic map construction
All the ltered SNPs/InDels were used to construct bin-map through MPR method. A total of 1354 bin markers covering 538,090 variants were identi ed on the 13 chromosomes ( Fig. 1). By using genotype data of the RIL population, we construct a high-density genetic map with a total genetic distance of 1295.45 cM.
The mapped bin per chromosome ranged from 58 (chr7) to 155 (chr3) with an average of 104.2 per chromosome (Fig. 2, Table 1). The density of bin markers in the whole genome was 0.98 cM/locus, covering an average physical length of 158.74 kb per bin. In addition to several gaps of more than 10 cM identi ed on chr2 (1), chr4 (1), chr5 (3), chr8 (1) and chr10 (2), the bin markers were distributed evenly along 13 chromosomes. The chromosome with the longest genetic length was chr2, which contained 120 bin markers covering a genetic length of 161.49 cM. Chr7 covered the shortest genetic length (55.01 cM) with 58 bin markers. A total of 1286 (95.0%) bin markers were less than 500 kb in length, and 40 bins covered physical length larger than 1 Mb. The largest bin located on chromosome 9 (c09b114), with a physical length larger than 5 Mb. (Addidtional le 2: Table S2).
To evaluate the quality of this high density genetic map, we investigated the colliearity between this genetic map and physical map. The dot plot of markers in the 13 linkage groups aligned well with the Zhongzhi No.13 reference chromosome, indicating excellent colliearity between genetic map and physical map (Fig.   3).

Phenotypic variation and correlation analysis
The phenotype data analysis showed that each trait varied among two parents and different RILs in both treatments ( Fig. 4; Table 2). The values of traits root length (RL), shoot length (SL) and fresh seedling weight (SW) of parent JHM were all lower than ZSB to some extent. The phenotypic distributions of mean showed continuous variations and transgressive segregations on both directions of the parents, suggesting polygenic inheritance of all traits in sesame. In the PEG osmotic condition, mean SL and SW of the RIL population were signi cant inhibited and reduced by 29.3% and 37.1% respectively when compared with the control condition, whereas the mean RL trait was slightly affected and reduced by 2.7%.
To better understand the responses of these traits to drought stress, we also investigated the relative phenotypic data of these three traits ( Fig. 5; Table 2). The mean values of relative SW (RSW) and SL (RSL) showed no differences between two parents and mean relative RL (RRL) of ZSB (0.88) was signi cantly higher than that of JHM (0.58), indicating that JHM was more sensitive to drought than ZSB at early seedling stage. In RILs, highly signi cant differences were also noted for all three relative parameters. Among them, RRL showed the highest heritability (75%) and the most remarkable variation, ranging part from 1 in two directions, indicating the inhibition or induction of osmotic stress on root growth in the population. This also suggested that the responses to drought stress in root of lines were mostly determined by their genotypes. The broad sense heritability of all traits under both conditions ranged from 34% to 88%, with the highest heritability values recorded for root length under PEG osmotic condition (88%), and relative seedling weight had the lowest value of 34% heritability ( Table 2). Correlations among seedling weight, shoot length and root length in both conditions were also surveyed (  The correlation coe cients from experiment 1 and experiment 2 are shown in the rst and second rows, respectively * and ** indicates signi cance at P = 0.05 and 0.01, respectively QTL mapping Composite interval mapping was used to identi ed chromosome regions associated with seedling weight, shoot length and root length under control and PEG stress conditions. In total, 34 QTLs were detected for all the traits under both conditions including the relative traits values, and 13 of them could be detected under two experiments. These stable QTLs were mapped to six chromosomes.

QTL identi cation under control condition
Under control condition, eleven QTLs were identi ed as being associated with these traits (Table 4). qSLC1, which could be identi ed in both experiments, in uenced shoot length, located at chromosome 1, had the highest LOD value, explaining an average of 16.26% of the phenotypic variation. Three chromosome regions showed associated with SW. Among them, qSWC12 was mapped to chromosome 12 and expressed stably across the trials, contributing to higher SW through JHM alleles. For shoot length, besides qSLC1, another two stable QTLs were detected on chromosome 5 and 12. One of them, qSLC12, was mapped to the interval (c12b062-c12b072) on chromosome overlapped with the interval of qSWC12 (c12b069-c12b071) on chromosome 12. qSLC12 had the second highest LOD value and explained up to 9.38% of phenotypic variation in SL. qSLC1 gained favourable allele from ZSB, while qSLC5 and qSLC12 were associated with increased SL through JHM alleles.
A total of four root length QTLs were detected under control condition. One stable QTL (qRLC4) located on chromosome 4 were detected in two trials (Table 4), had an average LOD value of 4.73 and explained at least 8.33% of phenotypic variation. The ZSB allele of qRLC4 was associated with longer root length.

QTL identi cation under PEG stress condition
Eleven chromosome regions were found to be associated with SW, SL and RL under PEG stress condition (SWP, SLP and RLP) ( Table 4). For SWP, two QTLs were identi ed and mapped on chromosome 3 and 9 (qSWP3 and qSWP9), but none of them could be detected across two trials. Four chromosome regions were detected to be associated with SL. Two of them located on chromosome 9. One QTL located on chromosome 8 (qSLP8) were steadily expressed in both trials and explained an average of 5.86% of phenotypic variation. qSLP8 increased SL under osmotic condition through the ZSB allele.
Of the four QTLs associated with RL under PEG stress condition, three could be identi ed in two experiments (Table 4). Among them, qRLP12 had the highest LOD value and strongest effect on RL, explaining up to 14.46% of phenotypic variation. The other two QTLs were located on chromosome 1 (qRLP1) and 7 (qRLP7), explaining at least 4.36% and 4.14% of phenotypic variation respectively. The ZSB alleles of these three stable QTLs were associated with longer root under PEG stress condition.
QTL identi cation for drought tolerance index DTI was de ned as relative traits value in this study. Twelve chromosome regions were identi ed to be associated with relative seedling weight (RSW), relative shoot length (RSL) and relative root length (RRL) ( Table 4; Additional le1: Figure S1). For RSL, qRSL1-2 detected in both trials had the strongest effect on RSL, explaining an average of 10.82% phenotypic variation, with LOD scores > 5. JHM contributed the positive alleles for qRSW5 and qRSL1-2. Four QTLs were identi ed for RSW in one of the trials. No stable QTL was detected for RSW.
Among the ve chromosome regions associated with RRL, qRRL12 had the strongest effect and explained, on average, 13.37% of the phenotypic variation in two trials, with an average LOD score of 7.81. This QTL was mapped to the same chromosome region as qRLP12 ( Fig. 6; Table 4). Besides qRRL12, two QTLs, qRRL1 and qRRL7, contributed at least 6.42% and 5.14% of the phenotypic variation in two trials, respectively. All of these three stable QTLs were responsible for the increase of RRL through the ZSB allele.
Putative candidate genes of major QTLs for drought tolerance Since the relative traits values, such as RSW, RSL and RRL, re ect the response to drought stress, we identi ed putative candidate genes of major stable QTLs for drought tolerance. Four genomic regions associated with ve major QTLs (qRSL1-2, qRRL1, qRRL7 and qRRL12/qRLP12) were obtained according to the corresponding anking markers of each QTL. A total of 465 genes were detected in these four regions, among them, 347 genes were functional annotated while 118 genes were identi ed as unknown protein, hypothetical proteins or repeat sequences. All the genes were used to identi ed homologs in Arabidopsis thaliana to help with gene functional annotation. 360 of them could detect A. thaliana homologs. Variations in the gene region between two parents from the re-sequencing data were also examined.
For RSL, the interval of major QTL qRSL1-2 contained 92 genes, of which, 74 were annotated as known protein-encoding genes. 73 genes could nd at least one homolog in A. thaliana genome, whereas 19 genes did not hit any A. thaliana genes. Four genes, SIN_1008253, SIN_1008230, SIN_1008220 and SIN_1008205, encoding GTP binding domain, FAD/NAD(P)-binding domain, RNA recognition motif domain and Zinc nger (A20-type) respectively, were homologous to Arabidopsis genes whose functional roles in response to drought stress had been validated through mutagenesis or transgenic studies [24][25][26][27]. Two of them, SIN_1008253 and SIN_1008220 contained nonsynonymous substitutions in the exons between two parents (Additional le 1: Table S3). SIN_1008253, homologous with AT3G57180, which encodes a chloroplast localized protein YL1/BPG2 that involved in seedling shoot response to salt stress [24]. Loss of YL1 function caused Na + accumulation and hypersensitive phenotype of shoot under salt stress [24]. SIN_1008230 harbored one SNP and a 17bp InDel in the promoter region (at the position of <-1kb from start codon).
Overexpressing seedlings of AT5G67030 (the homologs of SIN_1008230) displayed enhanced drought tolerance than the wild-type plants under osmotic stress [25].

Discussion
Sesame is one of the most important oil crops worldwide and provides sorts of speci c lignins which are very good for human health. Drought is a major stress effecting sesame growth at seedling stage. Efforts have been made for drought tolerance QTL and candidate genes identi cation in recent years [5,17].
However, the drought tolerance related traits are quite complex and controlled by multiple genes and environmental factors. Very few genetic studies were performed for drought tolerance at sesame seedling stage, especially for root related traits and no QTL was reported. In this study, 180 RILs were used to identify genomic regions associated with PEG-induced drought tolerance traits at early seedling stage in sesame.
Seedling weight, shoot length and root length were evaluated under both control and stress conditions.
Through whole genome re-sequencing, a high-density genetic map was constructed. 13 stable QTLs associated with drought tolerance traits were detected.
In recent decades, with the development of next-generation sequencing (NGS) technology and the release of whole genome sequence of sesame, numerous SNP markers have been used for QTL genetic mapping and GWAS analysis. Bin-maps constructed by high-density SNP markers using NGS have been widely used in sesame QTL mapping [22,40]. In the present study, we re-sequenced each RIL with an average read depth of 4.18× and construct a high-density genetic map with 1354 bins on all 13 chromosomes. The average genetic length was 0.978 cM per bin, which was similar with previous maps [22,40]. Comparing with other sesame genetic maps developed by restriction-site associated DNA sequencing (RAD-seq) or GBS approach, WGRS strategy could acquire much more variants. We obtained 466,911 high-quality SNPs and 72,981 InDels for bin markers assignment, which was more than 30 times of 13,679 SNPs identi ed by GBS in Zhang et al. [22] and 11,924 SNPs detected by RAD-seq in Zhang et al. [21]. The high-density variants could improve the resolution of QTL and also help perform ne mapping of QTL controlling various agronomic traits.
For drought tolerance evaluation, due to the unpredictability of rainfall and soil heterogeneity of eld experiment, PEG solution has been widely employed as an alternative to simulate water shortage condition, especially for the trials at germination and seedling stage [4,6,13,17]. In sesame, Li et al. [17]  of the phenotypic variation. These results implied that this interval on chromosome 12 contains a major QTL that associated with root length and drought tolerance. To our knowledge, this is the rst report of QTLs that associated with root length in sesame. Since the root is the tissue responsible for water uptake, deep root system can help plants absorbing water from deep soil layers to avoid drought stress. This interval in ZSB could be useful for drought tolerance breeding in sesame. In addition to qRLP12, the interval of other two stable QTLs controlling root length under osmotic stress condition were both closely linked with that of QTLs associated with relative root length (qRLP1 and qRRL1, qRLP7 and qRRL7) (Fig. 6). This strongly indicated close relationship between traits RLP and RRL, which may enable selection for the complex drought tolerance trait through an easily observable related trait, such as relative root length and root length under stress condition.
The traits of above ground part also play important role for avoiding drought. Long shoot is conducive to rapid seedling formation and allows deeper sowing. In this study, eight QTLs were identi ed for SL under both conditions, four of them (qSLC1, qSLC5, qSLC12 and qSLP8) stably expressed in two trials. For SW, only one major QTL qSWC12 for SW under control condition was detected in two trials. We could not detect any stable QTL for trait SW under osmotic stress condition and trait RSW, which may be due to their relative low heritabilities. This indicated that seedling weight may not be suitable as selection trait for drought tolerance in seedling stage of sesame. qSWC12 and qSLC12 was both located on similar region of chromosome 12. QTLs for shoot length and seedling weight co-location on chromosomes have also been report in soybean [42]. The JHM alleles for qSWC12 and qSLC12 were related to increased SW and SL under control condition, indicating the parent with lower phenotypic values may also contain favorable alleles for traits SW and SL.
QTL cluster in this region could be raised by both linkage and pleiotropic effects.
The positions of mapped QTLs in this study were also compared with previous drought tolerance related loci. Li et al. [17] detected nine and 15 QTLs for drought and salt tolerance indexes by utilizing germination rate and fresh seedling weight as phenotypes through GWAS at germination stage. None of the loci associated with fresh seedling weight were located in the QTL regions identi ed in our study, but a QTL qSGR6.1 related to germination rate index under salt stress shared the same location as qRRL12. Thabet et al. [43] detected several SNPs associated with RL_DTI and germination percentage_DTI on barley chromosome 3H and 4H were very close to each other at the corresponding regions. Whether the interval of qRRL12 were associated with germination rate in our population need further experiments. Dossa et al. [18] identi ed four genomic regions associated with ve traits including stem length and yield related traits at owering stage of sesame with GWAS. All of them was located at the positions different from those of QTLs in our study, suggesting various genetic mechanisms of drought tolerance may exist between the germination and owering stages.
Candidate genes located within the genomic region of four major QTLs for drought tolerance were also identi ed in this study. For trait RSL, we found three genes in qRSL1-2 region may be putative candidates according to the previous studies of their homologs in A. thaliana [24][25][26]. It is noteworthy that one of them, SIN_1008253, which contained seven nonsynonymous SNPs in the exons between two parents, is homologous with a well-studied A. thaliana gene YL1 (AT3G57180) that in uenced seedling shoot growth under salt stress [24]. YL1 encodes a YqeH-type GTPase that is conserved in plant chloroplast and regulate the shoot response to salt stress through ABI4 [24]. SIN_1008220, encoding a peptide containing RNA recognition motif domain. Its homologs in A. thaliana (PSRP2, AT3G52150) is a chloroplast-localized ribosomal protein that play a role as a negative regulator on seedling growth under salt stress condition [26].
Another gene SIN_1008230, contained variations in its promter region between parents, is homologous with a zeaxanthin epoxidase gene in A. thaliana (AtZEP, AT5G67030) that is an important enzyme in ABA biosynthesis. ZEP-overexpressing plants showed more developed rosette leaves and lateral roots than wildtype plants under mannitol stress condition [25].
For trait RRL, ve putative candidate genes located on three major QTLs regions (qRRL1, qRRL7 and qRRL12/qRLP12) were identi ed in this study. SIN_1022428 was the only one gene which harbored variations that caused peptide change among all 55 annotated genes in qRRL12/qRLP12 region. The A. thaliana homolog (AT1G06490) of SIN_1022428 encodes callose synthase 7, which evolves in callose deposition in developing sieve elements during phloem formation. The 5-day-old seedling of AT1G06490 mutant developed much shorter roots than wild-type [38]. Two genes, SIN_1017975 and SIN_1017978, located at the interval of qRRL1, were annotated as WRKY transcription factor and AP2/ERF protein, whose homologs in A. thaliana were closely related to abiotic stress responses [29]. SIN_1017978, also named as AP2si16, was reported to be highly induced by drought stress condition in previous study [44]. In the region of major QTL qRRL7, SIN_1006102 and SIN_1011692 which encode eIF-like protein and multi antimicrobial extrusion protein were the two genes included synonymous SNP. SIN_1006102 was homologous with AT5G57870 that was proved to be involved in the responses to dehydration, salinity, and heat stress through mutagenesis analysis [32], while SIN_1011692-encoding protein shared high similarity with a detoxi cation e ux carriers (DTX)/multidrug and toxic compound extrusion (MATE) transporter AT5G65380 (DTX27).
Overexpression of a cotton homolog of DTX27 could enhance drought, salt, and cold stress tolerance in transgenic Arabidopsis plants [39]. The con rmation of candidate genes needs further ne mapping of these QTLs and functional validation experiments.

Conclusions
In conclusion, we performed rst QTL mapping of drought tolerance related traits using a RIL population and identi ed root length QTLs for the rst time in sesame. By using PEG treatment, inheritances of three traits including seedling weight, shoot length and root length were interpreted. Root length had the largest broad sense heritability, while seedling weight occupied the lowest position. Using WGRS technique, we identi ed a total of 13 stable QTLs for drought tolerance related trait at seedling stage of sesame. Four of them explained more than 10% phenotypic variation and had an LOD score larger than 6. The current study was the rst report for QTL identi cation for root length, which is considered as one of the most important traits for drought tolerance in plant. A genomic interval on chromosome 12 contained two QTLs associated to RLP and RRL. We identi ed seven candidate genes that may be related to drought tolerance according to their gene descriptions and variations between parents. The major QTL regions and linked markers can provide potential genetic resources for molecular marker-assisted selection and further cloning of functional genes for drought tolerance in sesame.

Plant materials and phenotyping
The Sesamum indicum L. cultivars Jinhuangma (JHM) and Zhushanbai (ZSB) were landraces collected at Jiangxi province and Hubei province respectively of China. The seeds of JHM and ZSB used in this study were obtained from sesame germplasm reservoir of Nanchang Branch of National Center of Oilcrops Improvement, Jiangxi Academy of Agricultural Sciences, China. A set of 180 F 9 recombinant inbred lines developed by single seed descent from the cross of JHM and ZSB was used for QTL mapping. All the plants were grown in a nylon net house to prevent cross-pollination caused by insects.
Polyethylene glycol (PEG)-simulated drought stress trials were performed by using both parent and RIL population in two years. For each line, 50 mature seeds were germinated on two layers of lter paper in each plastic container (

Linkage map construction
The SNPs or InDels were ltered according to three criteria: (1) missing data rate < 20%; (2) minor allele frequency (MAF) > 0.2; and (3) loci that were homozygous in both parents and heterozygous in less than 15% of the RILs. Filtered SNPs/InDels were used to construct bin-map through maximum parsimonious inference of recombination (MPR) method described by Xie et al. [48]. First, all variants were re-ltered by permutations involving resampling of SNPs/InDels windows and then inferred by Bayesian method. The genotype of each locus in RILs was determined assisted by hidden Markov model. Consecutive SNPs/InDels sites with the same genotype as one parent were assembled into a block. A recombination event was de ned as a transition between two blocks with different genotypes. The R/qtl package (http://rqtl.org/) was used to construct the linkage map.

QTL analysis
The broad sense heritability was estimated with the formula: , where σ 2 g is the estimated genetic variance and σ 2 e is experimental error. QTLs for each trait in each experiment were identi ed by composite interval mapping (CIM) method using Window QTL Cartographer v2.5 [49] (http://statgen.ncsu.edu/qtlcart/WQTLCart.htm). The standard model 6 and a window size of 10 cM was applied. The level of signi cance was determined with 1000 permutations, with a con dence level of 95%.
The LOD score for declaring a QTL was 2.5 or above. MapChart software was used to construct the graphical representation of QTL positions [50].

Declarations
Acknowledgements Not applicable.

Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. Figure 1 Genotype data of individuals of the RIL population, red and blue represent Zhushanbai and Jinhuangma, green represents heterzygote.