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 filtering, a total of 466,911 high-quality SNPs and 72,981 InDels (Insertion and deletions) were identified 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 file 1: Table S1).
High density genetic map construction
All the filtered SNPs/InDels were used to construct bin-map through MPR method. A total of 1354 bin markers covering 538,090 variants were identified 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.
Table 1 Characteristics of the high-density genetic map
Chr.a
|
Length(cM)
|
No. markers
|
No. bins
|
Average bin interval (cM)
|
Max interval (cM)
|
chr1
|
89.539
|
61646
|
125
|
0.722
|
4.528
|
chr2
|
161.486
|
38457
|
120
|
1.357
|
11.684
|
chr3
|
128.823
|
76256
|
155
|
0.837
|
5.969
|
chr4
|
108.981
|
26753
|
77
|
1.434
|
25.697
|
chr5
|
116.469
|
43314
|
97
|
1.213
|
15.049
|
chr6
|
114.018
|
60058
|
136
|
0.845
|
6.339
|
chr7
|
55.011
|
12058
|
58
|
0.965
|
3.833
|
chr8
|
111.070
|
45042
|
119
|
0.941
|
10.355
|
chr9
|
91.418
|
37416
|
114
|
0.809
|
9.503
|
chr10
|
123.251
|
36959
|
105
|
1.185
|
29.912
|
chr11
|
65.639
|
40345
|
85
|
0.781
|
6.339
|
chr12
|
71.325
|
31949
|
77
|
0.938
|
9.503
|
chr13
|
58.423
|
27837
|
86
|
0.687
|
3.491
|
Whole
|
1295.453
|
538090
|
1354
|
0.978
|
10.939
|
aChr., Chromosome
In addition to several gaps of more than 10 cM identified 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 file 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 significant 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 significantly higher than that of JHM (0.58), indicating that JHM was more sensitive to drought than ZSB at early seedling stage. In RILs, highly significant 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).
Table 2 Phenotypic data and heritability of six traits in parents and RIL population
Traits
|
Experiment
|
|
ZSB
|
JHM
|
RIL mean
|
RIL range
|
h2
|
Seedling weight (mg)
|
1
|
Control
|
73.35
|
62.10
|
60.41
|
46.56-76.50
|
0.62
|
|
2
|
Control
|
73.13
|
62.75
|
66.36
|
39.40-79.88
|
|
|
1
|
PEG
|
40.65
|
35.20
|
38.87
|
29.67-58.08
|
0.55
|
|
2
|
PEG
|
41.95
|
35.25
|
40.89
|
12.50-60.93
|
|
Shoot length (cm)
|
1
|
Control
|
5.04
|
4.88
|
4.61
|
3.77-5.47
|
0.73
|
|
2
|
Control
|
5.22
|
4.94
|
4.84
|
3.93-6.56
|
|
|
1
|
PEG
|
3.18
|
3.16
|
3.27
|
2.67-3.86
|
0.71
|
|
2
|
PEG
|
3.19
|
3.09
|
3.41
|
2.55-3.98
|
|
Root length (cm)
|
1
|
Control
|
7.48
|
5.78
|
6.63
|
4.94-8.18
|
0.71
|
|
2
|
Control
|
7.46
|
5.86
|
6.15
|
3.91-7.43
|
|
|
1
|
PEG
|
6.48
|
3.30
|
6.34
|
3.24-8.67
|
0.88
|
|
2
|
PEG
|
6.59
|
3.49
|
6.09
|
3.19-8.75
|
|
Relative seedling weight
|
1
|
|
0.55
|
0.57
|
0.64
|
0.53-0.95
|
0.34
|
|
2
|
|
0.57
|
0.56
|
0.65
|
0.16-1.02
|
|
Relative shoot length
|
1
|
|
0.63
|
0.65
|
0.71
|
0.56-0.87
|
0.70
|
|
2
|
|
0.61
|
0.62
|
0.71
|
0.53-0.91
|
|
Relative root length
|
1
|
|
0.87
|
0.57
|
0.96
|
0.61-1.33
|
0.75
|
|
2
|
|
0.88
|
0.59
|
0.97
|
0.58-1.41
|
|
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 significant relationships with RLP (significant at P = 0.05) in either of two trials. All the identical traits showed significant correlation between control and stress condition, with the average correlation coefficients ranging from 0.347 for SL to 0.472 for RL.
Table 3 Correlation coefficents for six traits in the RIL population
|
SW
|
SL
|
RL
|
SWP
|
SLP
|
RLP
|
SW
|
1
|
|
|
|
|
|
SL
|
0.193**
0.390**
|
1
|
|
|
|
|
RL
|
0.295**
0.381**
|
0.249**
0.475**
|
1
|
|
|
|
SWP
|
0.776**
0.461**
|
-0.008
0.042
|
0.240**
0.086
|
1
|
|
|
SLP
|
0.031
0.288**
|
0.303**
0.390**
|
0.102
0.218**
|
0.091
0.356**
|
1
|
|
RLP
|
0.001
0.117
|
0.172*
0.380**
|
0.454**
0.489**
|
-0.099
0.072
|
0.414**
0.509**
|
1
|
The correlation coefficients from experiment 1 and experiment 2 are shown in the first and second rows, respectively
* and ** indicates significance at P = 0.05 and 0.01, respectively
QTL mapping
Composite interval mapping was used to identified 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 identification under control condition
Under control condition, eleven QTLs were identified as being associated with these traits (Table 4). qSLC1, which could be identified in both experiments, influenced 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,
Table 4 QTLs identified for drought tolerance traits
Traits
|
Chromosome
location
|
Locus
|
Flanking Markers
|
Experiment 1 (2018)
|
|
Experiment 2 (2019)
|
LODa
|
AEb
|
R2(%)
|
|
LODa
|
AEb
|
R2(%)
|
SW-control
|
Chr2
|
qSWC2
|
c02b067-c02b073
|
|
|
|
|
3.22
|
1.71
|
6.54
|
|
Chr8
|
qSWC8
|
c08b086-c08b090
|
3.75
|
1.62
|
7.17
|
|
|
|
|
|
Chr12
|
qSWC12
|
c12b069-c12b071
|
5.56
|
1.99
|
10.82
|
|
2.98
|
1.56
|
6.05
|
SL-control
|
Chr1
|
qSLC1
|
c01b092-c01b100
|
9.06
|
-0.14
|
16.38
|
|
8.72
|
-0.17
|
16.14
|
|
Chr5
|
qSLC5
|
c05b087-c05b094
|
3.05
|
0.08
|
5.02
|
|
2.87
|
0.10
|
4.87
|
|
Chr8
|
qSLC8
|
c08b065-c08b070
|
3.54
|
0.08
|
5.86
|
|
|
|
|
|
Chr12
|
qSLC12
|
c12b062-c12b072
|
5.34
|
0.10
|
9.06
|
|
5.30
|
0.13
|
9.38
|
RL-control
|
Chr1
|
qRLC1
|
c01b086-c01b087
|
|
|
|
|
5.23
|
-0.20
|
10.46
|
|
Chr4
|
qRLC4
|
c04b040-c04b046
|
4.91
|
-0.21
|
9.12
|
|
4.54
|
-0.19
|
8.33
|
|
Chr6
|
qRLC6
|
c06b129-c06b135
|
|
|
|
|
4.57
|
0.18
|
8.46
|
|
Chr10
|
qRLC10
|
c10b076-c10b080
|
4.62
|
-0.20
|
8.40
|
|
|
|
|
SW-PEG
|
Chr1
|
qSWP1
|
c01b003-c01b010
|
4.20
|
-1.23
|
8.20
|
|
|
|
|
|
Chr3
|
qSWP3
|
c03b116-c03b119
|
|
|
|
|
3.16
|
1.39
|
6.61
|
|
Chr9
|
qSWP9
|
c09b031-c09b040
|
|
|
|
|
3.57
|
-1.62
|
7.50
|
SL-PEG
|
Chr1
|
qSLP1
|
c01b032-c01b035
|
4.27
|
-0.07
|
8.26
|
|
|
|
|
|
Chr8
|
qSLP8
|
c08b055-c08b063
|
3.06
|
0.06
|
5.88
|
|
3.04
|
0.07
|
5.83
|
|
Chr9
|
qSLP9-1
|
c09b015-c09b021
|
|
|
|
|
3.76
|
0.10
|
7.26
|
|
Chr9
|
qSLP9-2
|
c09b031-c09b033
|
|
|
|
|
5.56
|
-0.12
|
11.00
|
RL-PEG
|
Chr1
|
qRLP1
|
c01b062-c01b070
|
4.97
|
-0.31
|
7.97
|
|
4.36
|
-0.31
|
6.77
|
|
Chr6
|
qRLP6
|
c06b054-c06b060
|
3.28
|
0.26
|
5.15
|
|
|
|
|
|
Chr7
|
qRLP7
|
c07b030-c07b036
|
4.37
|
-0.29
|
6.61
|
|
4.14
|
-0.30
|
6.40
|
|
Chr12
|
qRLP12
|
c12b032-c12b036
|
7.17
|
-0.38
|
11.85
|
|
8.78
|
-0.45
|
14.46
|
RSW
|
Chr5
|
qRSW5-1
|
c05b019-c05b029
|
3.74
|
0.01
|
6.80
|
|
|
|
|
|
Chr5
|
qRSW5-2
|
c05b071-c05b074
|
|
|
|
|
3.01
|
0.02
|
6.33
|
|
Chr6
|
qRSW6
|
c06b041-c06b045
|
4.51
|
0.01
|
8.29
|
|
|
|
|
|
Chr12
|
qRSW12
|
c12b061-c12b072
|
|
|
|
|
3.56
|
-0.02
|
7.47
|
RSL
|
Chr1
|
qRSL1-1
|
c01b035-c01b049
|
|
|
|
|
3.56
|
-0.17
|
6.39
|
|
Chr1
|
qRSL1-2
|
c01b109-c01b113
|
5.37
|
0.02
|
9.48
|
|
6.67
|
0.02
|
12.15
|
|
Chr11
|
qRSL11
|
c11b044-c11b051
|
3.74
|
0.02
|
6.36
|
|
|
|
|
RRL
|
Chr1
|
qRRL1
|
c01b052-c01b063
|
4.12
|
-0.04
|
6.56
|
|
3.92
|
-0.04
|
6.42
|
|
Chr3
|
qRRL3-1
|
c03b043-c03b055
|
|
|
|
|
3.65
|
-0.04
|
5.97
|
|
Chr3
|
qRRL3-2
|
c03b102-c03b113
|
4.29
|
-0.04
|
6.82
|
|
|
|
|
|
Chr7
|
qRRL7
|
c07b020-c07b028
|
4.45
|
-0.04
|
7.09
|
|
3.19
|
-0.04
|
5.14
|
|
Chr12
|
qRRL12
|
c12b032-c12b036
|
6.26
|
-0.05
|
10.26
|
|
9.36
|
-0.07
|
16.47
|
aLOD: likelihood of the odds
bAdditive effect: positive and negtive indicated ZSB and JHM allele produced larger value respectively.
cE followed by 1 or 2 designate two independent experiments
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 identification 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 identified 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 identified 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 identification for drought tolerance index
DTI was defined as relative traits value in this study. Twelve chromosome regions were identified to be associated with relative seedling weight (RSW), relative shoot length (RSL) and relative root length (RRL) (Table 4; Additional file1: 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 identified for RSW in one of the trials. No stable QTL was detected for RSW.
Among the five 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, reflect the response to drought stress, we identified putative candidate genes of major stable QTLs for drought tolerance. Four genomic regions associated with five major QTLs (qRSL1-2, qRRL1, qRRL7 and qRRL12/qRLP12) were obtained according to the corresponding flanking 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 identified as unknown protein, hypothetical proteins or repeat sequences. All the genes were used to identified 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 find 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 finger (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-27]. Two of them, SIN_1008253 and SIN_1008220 contained nonsynonymous substitutions in the exons between two parents (Additional file 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].
For RRL, the physical regions of three major QTLs (qRRL1, qRRL7 and qRRL12/qRLP12) on chromosome 1,7 and 12 encompassed 157,161 and 55 predicted genes respectively. 109, 95 and 38 of them were described as functional genes. A total of 287 genes could identify homologous genes in A. thaliana. 12 genes of them, encoding cellulose synthase, eukaryotic translation initiation factor (eIF)-like protein, subtilisin-like protease, leucine aminopeptidase/peptidase B, ELO family member, multi antimicrobial extrusion protein, aquaporin-like protein, aquaporin-like protein, WRKY transcription factor, AP2/ERF protein and ATP-dependent RNA helicase DEAD-box, had been reported to regulate responses to drought stress in Arabidopsis [28-39]. Four of these genes comprised nonsynonymous SNPs in the coding region between two parents, including SIN_1017978 encoding AP2/ERF domain in qRRL1 region, SIN_1006102 and SIN_1011692 encoding eIF-like protein and multi antimicrobial extrusion protein in qRRL7 region and SIN_1022428 encoding callose synthase 7 in qRRL12/qRLP12 region (Additional file 3: Table S3). Two genes, SIN_1022428 and SIN_1017975 (within qRRL1 region), contained variations in their promoter regions.