Construction of a High Density Genetic Map and QTL Mapping of PEG-Induced Drought Tolerance at Seedling Stage in Sesame Using Whole Genome Re-Sequencing

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, very 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 identied under both non-stress and water stress conditions. were evaluated under and PEG-induced osmotic conditions at stage in a from cross of Zhushanbai and Jinhuangma Signicant variation and high broad sense heritability were observed for all traits except SW under stress condition in the population. With this population, a high-density 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 three traits under both conditions and their relative values, and 13 stable QTLs associated with seven traits could be 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. Root length related QTLs were rst identied in One region on chromosome 12 contained two major QTLs related to RL under osmotic condition and relative RL. for high resolution mapping of QTLs in sesame [18–22]. QTL identication 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 identied three genes that may be related to the regulation of root length under osmotic stress. 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.


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 [5,6,12].
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 [13][14][15]. Mensah et al. [13] 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 few studies have been conducted yet on genetic analysis of drought tolerance. Li et al. [16] 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 [17]. 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 [17]. 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 resequencing. These methods have been successfully utilized for high resolution mapping of QTLs in sesame [18][19][20][21][22].
To further explore the genetic foundation of sesame drought tolerance, in the present study, we developed a highdensity 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 [19] 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. Chromosome (chr) 3 harbored the largest number of SNPs and InDels, while chromosome 7 contained the fewest number of variants. The density of the SNP and InDel loci in the genome was 1788.15/Mb and 280.32/Mb respectively (Additional le 1: 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). In RILs, highly signi cant differences were also noted for all three relative parameters. Among them, RRL showed remarkable variation, ranging part from 1 in two directions, indicating the inhibition or induction of osmotic stress on root growth in the population. 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 ( Table 3). For control condition, all the three traits (herein abbr. SWC, SLC and RLC, repectively) were positively correlated with each other in the two trials. Correlations among SW, SL and RL under PEG stress condition (herein abbr. SWP, SLP and RLP, respectively) were weaker than that under control condition. RLP were positively correlated with SLP in one or two experiments, but SWP had no signi cant relationships with RLP (signi cant at P = 0.05) in either of two trials. All the identical traits showed signi cant correlation between control and stress condition, with the average correlation coe cients ranging from 0.347 for SL to 0.472 for RL. 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, 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
Drought tolerance index 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).
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.

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,16]. 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 [21,23].
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 [21,23]. Comparing with other sesame genetic maps developed by restriction-site associated DNA sequencing (RAD-seq) or GBS approach, WGS 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. [23] and 11,924 SNPs detected by RAD-seq in Zhang et al. [21]. The highdensity 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 seedling stage [6,16,24]. In sesame, Li et al. [16] used relative values of germination rate and fresh weight (the ratio of the traits value under stress conditions to the same traits under stress free conditions) to evaluate the drought tolerance of sesame at early seedling stage. For drought tolerance at sesame owering stage, Dossa et al. [17] investigated wilting level of the whole plant, stem length, and some yield related traits. However, no comprehensive genetic study has been conducted yet for the root traits under desiccation stress condition in sesame. In this study, three characters including SW, SL and RL were measured in the RIL population under control and PEG stress condition.
All the measured traits showed signi cant genetic variation under both conditions, which enable the identi cation of QTLs associated with traits in sesame. We found high level of broad heritability of the trait RL under stress condition, indicating the selection of this trait could be effective in the breeding program. Relative performances of trait values were usually used as indicators for the evaluation of stress tolerance. We also investigated SW, SL and RL in current study and the results showed that relative RL distributed part from 1 in both directions, suggesting the PEG solution treatment in our study only inhibited the root growth of some lines. This also indicated that the different performances of root length in the RILs depended on genotypic variation.
In this study, QTLs were identi ed for all the traits under control and osmotic stress conditions. For root length, we found that in ZSB, both qRLP12 and qRRL12 were mapped to the same interval on chromosome 12, associated with increased RL under PEG stress condition and RRL. qRLP12 and qRRL12 had the highest LOD values of all the QTLs associated with RLP and RRL respectively, and explained up to 14.46 and 16.67% 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 root length.
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 [25]. 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.
Based on the QTL mapping results, we examined the genes located within the genomic region of qRLP12/qRRL12, which had the largest effect on drought tolerance in this study. This region was related with both RL under PEG stress condition and relative RL. 55 annotated genes were found between the boundary markers of the interval of qRLP12/qRRL12 in the sesame reference genome (Additional le 3: Table S3). Then, we identi ed homologs of these genes in Arabidopsis thaliana through BLASTX using the coding sequences of these 55 genes. Of them, 43 could nd homologs in A. thaliana, and three genes may be identi ed as candidates. SIN_1022439 encodes an uncharacterized peptide, which is homologous with AT1G17120 (RASD1, Responsiveness to ABA salt and drought 1). It had been reported that RASD1 was expressed predominantly in the vascular system, and rasd1 mutant lines were veri ed to be drought intolerant and less sensitive to exogenous ABA and NaCl during germination and root growth than wild-type plants [26]. SIN_1022428, homologous with AT1G06490, encoding callose synthase 7 (CalS7), 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 [27]. SIN_1022433 is another candidate encoding subtilisin-like protease. Its homolog in A. thaliana (AT2G04160, AIR3) was reported to be involved in lateral root morphogenesis and induced by NAC genes which regulated multiple abiotic stress responses [28,29]. The con rmation of candidate genes needs further ne mapping of this QTL.

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 three genes that may be related to the regulation of root length under osmotic stress. 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 crosspollination 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. [33]. 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 [34] (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 [35].        Bars in each chromosome represent bin-markers. Vertical bars of each QTL represent con dence intervals of 1-LOD and are located compared to the marker position on each chromosome.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.