Next-generation sequencing identified genomic region and diagnostic markers for resistance to bacterial wilt on chromosome B02 in peanut (Arachis hypogaea L.).

Bacterial wilt, caused by Ralstonia solanacearum, is a devastating disease affecting over 350 plant species. A few peanut cultivars were found to possess stable and durable bacterial wilt resistance (BWR). Genomics-assisted breeding can accelerate the process of developing resistant cultivars by using diagnostic markers. Here, we deployed sequencing-based trait mapping approach, QTL-seq, to discover genomic regions, candidate genes and diagnostic markers for BWR in a recombination inbred line population (195 progenies) of peanut. The QTL-seq analysis identified one candidate genomic region on chromosome B02 significantly associated with BWR. Mapping of newly developed single nucleotide polymorphism (SNP) markers narrowed down the region to 2.07 Mb and confirmed its major effects and stable expressions across three environments. This candidate genomic region had 49 nonsynonymous SNPs affecting 19 putative candidate genes including seven putative resistance genes (R-genes). Two diagnostic markers were successfully validated in diverse breeding lines and cultivars and could be deployed in genomics-assisted breeding of varieties with enhanced BWR.

Introduction

The soil‐borne bacterium Ralstonia solanacearum causes one of the most devastating diseases, bacterial wilt, over 350 plant species, including major crops such as tobacco, tomato, potato, pepper, eggplant, banana and peanut (Salanoubat et al., 2002). R. solanacearum is found worldwide and could survive many years in soils (Hayward, 1991). This pathogen usually infects plants through roots and then spreads through the vascular system to aerial parts. Plants will show wilting symptoms and die if these bacteria multiply to high level (Genin, 2010). It is difficult to eradicate R. solanacearum because of its wide host range and long‐term survival ability.

During the plant–pathogen co‐evolution, plants developed quantitative resistance (French et al., 2016; Yang et al., 2017) to bacterial wilt disease. The quantitative resistance could inhibit the multiplication of R. solanacearum in plants, which usually leads to a reduction in symptom but not the absence of disease. So far, only two resistance genes (R‐genes) for bacterial wilt resistance (BWR) have been cloned from model plant Arabidopsis thaliana. One major quantitative trait locus (QTL) was identified and cloned in A. thaliana, leading to identification of the ERECTA gene which encodes an LRR receptor‐like kinase (RLK) whose up‐regulated expression limited the pathogen growth in host thus avoiding wilting symptoms (Godiard et al., 2003). The other gene for BWR cloned in A. thaliana is RRS1‐R, encoding a TIR‐NBS‐LRR (TNL) class disease resistance protein and its expression could also limit the pathogen population in host (Deslandes et al., 2002). Further studies suggested that the RRS1‐R protein physically associates with RPS4 protein to trigger host resistance (Williams et al., 2014). Interestingly, when RRS1‐R/RPS4 was transferred into tomato, transgenic lines exhibited resistance to R. solanacearum (Narusaka et al., 2013). In addition, Zhang et al. (2017) reported that the overexpression of the peanut NBS‐LRR gene AhRRS5 in tobacco enhanced BWR of transgenic lines. These results indicated the potential use of R‐genes to control bacterial wilt across species.

The cultivated peanut (Arachis hypogaea L.), also known as groundnut, is an important legume crop grown in more than 100 countries and consumed worldwide in form of high‐quality cooking oil, edible nut, peanut butter and candy, etc. (Pandey et al., 2012; Varshney et al., 2013). In 2016, the global production of peanut was more than 43.98 million tonnes (FAOSTAT). China and India are the top two producers with production of 16.62 and 6.86 million tonnes in 2016, respectively. Bacterial wilt is a devastating disease in peanut production of many countries such as China, Vietnam, the Philippines, Indonesia, Malaysia, Thailand, India, Sri Lanka and Uganda. For example, bacterial wilt has been found in most of the 13 main peanut producing provinces in China, and it could cause up to 50%–100% yield losses (Jiang et al., 2017). The international Groundnut Bacterial Wilt Working Group (GBWWG) workshop was organized in 1998 to develop effective strategy for controlling losses from bacterial wilt disease. Deployment of genetic resistance by breeding peanut cultivars with stable resistance was found to be the most effective approach for controlling bacterial wilt disease. Dozens of resistant cultivars were developed through conventional breeding and deployed widely for peanut production (Mallikarjuna and Varshney, 2014). Their resistances were found to be the most stable and durable among plant species.

The information on genetic behaviour and trait mapping of BWR in peanut is still lacking and hence needs more efforts for understanding trait genetics, gene discovery and developing diagnostic markers. Availability of diagnostic markers would be beneficial to transfer BWR into more elite cultivars using genomics‐assisted breeding (GAB) in less time and more precisely. Two major QTLs, qBW‐1 and qBW‐2, for BWR in an F2 population (Yueyou 92 × Xinhuixiaoli) and only qBW‐1 could be confirmed in its recombination inbred line (RIL) population (Zhao et al., 2016). Due to the similarity between A and B subgenome of the cultivated peanut (AABB, 2n = 4x = 40), the big genome size (~2.7 GB) and the low polymorphism of molecular markers, it was challenge, time‐consuming and labour‐intensive to cover the whole genome with traditional molecular markers. For example, Zhao et al. (2016) screened 1343 SSR markers and 180 Kompetitive Allele Specific PCR (KASP) markers; however, the constructed linkage map had just 237 marker loci that were obviously not enough to cover the peanut genome. This study provided limited information of genomic locations of these QTLs without reporting any diagnostic marker for BWR. Therefore, a more effective strategy is needed for rapidly mapping QTLs for BWR in the other varieties, which could be further pyramided through genomics‐assisted breeding (GAB) to accelerate the process of breeding elite varieties with enhanced BWR.

Based on bulked segregant analysis (BSA; Michelmore et al., 1991) and next‐generation sequencing, the QTL‐seq approach (Takagi et al., 2013) was developed to rapidly identify genomic regions for trait of interest in plants. This approach required resequencing of three samples, that is two bulks segregating for trait of interest and one parent of the mapping population, and is therefore cost‐effective and more precise (Takagi et al., 2013). This approach has been proven very successful in identification of genomic regions for resistance to Fusarium wilt and sterility mosaic disease in pigeonpea (Singh et al., 2016), resistance to rust and late leaf spot in peanut (Pandey et al., 2017) and resistance to ascochyta blight in chickpea (Deokar et al., 2018). The present study reports deployment of sequencing‐based trait mapping approach, QTL‐seq, that helped in successful discovery of genomic regions, candidate genes and diagnostic markers for BWR in a RIL population developed from the cross of Yuanza 9102 × Xuzhou 68‐4.

Results

Phenotypic variations and construction of extreme bulks for bacterial wilt resistance

The percentage of plants which were not killed by the bacterial wilt disease until preharvest stage was used to assess the resistance, referred as survival rate hereafter. The survival rate of the resistant parent Yuanza 9102 (RP) was significantly higher than that of the susceptible parent Xuzhou 68‐4 (SP) (Figure 1a–c). The survival rates of the 195 RILs showed continuous distribution with two peaks across three environments (Figure 1b), indicating the existing of major QTLs for BWR. Based on the mean values of survival rates in the three environments, 15 RILs with lowest survival rates (2.15%–11.54%) and 15 RILs with highest survival rates (92.96%–100%) were selected to prepare the susceptible bulk (SB) and resistant bulk (RB), respectively (Figure 1c).

Figure 1

Phenotypic variations and construction of the extreme bulks for bacterial wilt resistance. (a) The difference of survival rate of susceptible parent Xuzhou 68‐4 and resistant parent Yuanza 9102 in the field disease nursery. (b) Frequency distribution for survival rates in RIL population. The y‐axis represented density, while the x‐axis represented survival rates. The dotted lines represented the survival rates of Xuzhou 68‐4, and the dashed lines represented the survival rates of Yuanza 9102. (c) Phenotypic variations among the RILs selected for developing extreme bulks for bacterial wilt resistance. Based on mean values of three environments, 15 RILs with lowest survival rates and 15 RILs with highest survival rates were used to construct susceptible and resistant bulks.

Whole‐genome resequencing, mapping of reads and identification of SNPs

Whole‐genome resequencing data, including 544.70 million reads (82.25 Gb) for SP, 456.89 million reads (68.99 Gb) for RP, 739.71 million reads (111.70 Gb) for SB and 729.63 million reads (110.17 Gb) for RB (Table S1), were generated using Illumina NovaSeq platform and analysed using the pipeline shown in Figure 2. Using the SP assembly as reference, the reads of the two extreme bulks were mapped and achieved 93.04% coverage and 21.41 X read depth for SB while 93.15% coverage and 20.96 X read depth for RB (Table S1). The comprehensive sequence analysis between the susceptible and resistant bulks identified 164 522 genome‐wide SNPs (Table S2). On the other hand, the reads of the two bulks were also mapped to the RP assembly and achieved 92.98% coverage and 20.40 X read depth for SB while 93.09% coverage and 20.96 X read depth for RB (Table S1). A total of 243 380 genome‐wide SNPs were identified between the extreme bulks (Table S2). The intersection of the two sets of SNPs derived from two different reference assemblies was 75 080.

Figure 2

QTL‐seq approach used for mapping bacterial wilt resistance in peanut.

Candidate genomic regions for bacterial wilt resistance

For each genome‐wide SNP identified above, the SNP‐index was firstly calculated for SB and RB (Figure S1). By subtracting SNP‐index of RB from SNP‐index of SB, the ΔSNP‐index were calculated for each genome‐wide SNPs. Based on the sliding window analysis for SNP‐index and ΔSNP‐index plots, genomic regions with SNP‐index significantly deviated from 0.5 and ΔSNP‐index significantly deviated from 0 (zero) were identified as candidates controlling BWR.

With the SP assembly as reference, a major peak on chromosome B02 was identified for BWR at a statistical confidence of P < 0.01 (Figure 3e–h, Figures [Link], [Link], [Link]), spanning a 3.35 Mb (3.25–6.60 Mb) interval. The ΔSNP‐index was negative, indicating that more alleles were from nonreference parent Yuanza 9102 in RB (Table 1). However, with the RP assembly as reference, a 2.85 Mb (13.55–16.40 Mb) interval on chromosome B09 was also identified for BWR, in addition to a similar (3.25–6.90 Mb) interval on chromosome B02 (Figure 3a–d, Figures [Link], [Link], [Link]). The ΔSNP‐index of the two genomic regions was positive, indicating that more alleles were from reference parent Yuanza 9102 in RB (Table 1).

Figure 3

Genome‐wide overview of the genomic regions identified for bacterial wilt resistance. (a) Pseudomolecules of reference genome Arachis duranensis and A. ipaensis (b) Genome‐wide density of annotated NBS‐LRR genes. (c) ΔSNP‐index plot using the Yuanza 9102 assembly as reference, from outside to inside: upper probability values at 99% confidence, upper probability values at 95% confidence, ΔSNP‐index, lower probability values at 95% confidence, lower probability values at 99% confidence. (d) SNP‐index plots of susceptible bulk (green) and resistant bulk (red) using the Yuanza 9102 assembly as reference. (e) ΔSNP‐index plot using the Xuzhou 68‐4 assembly as reference, from outside to inside: upper probability values at 99% confidence, upper probability values at 95% confidence, ΔSNP‐index, lower probability values at 95% confidence, lower probability values at 99% confidence. (f) SNP‐index plots of susceptible bulk (green) and resistant bulk (red) using the Xuzhou 68‐4 assembly as reference. (g) Physical position of QTL identified using SSR‐based genetic linkage map. (h) The highlights of the combined genomic regions listed in Table 1.

Table 1

Genomic regions identified for bacterial wilt resistance

Reference assembly	Chr§	Genomic region (Mb)	Length (Mb)	ΔSNP‐index¶	U99††	L99‡‡	Allele source
SP*	B02	3.25–6.60	3.35	−0.87	0.63	−0.62	Yuanza 9102
RP†	B02	3.25–6.90	3.65	0.87	0.63	−0.63	Yuanza 9102
	B09	13.55–16.40	2.85	0.66	0.63	−0.63	Yuanza 9102
Combined‡	B02	3.25–6.60	3.35				Yuanza 9102
	B09	13.55–16.40	2.85				Yuanza 9102

The susceptible parent Xuzhou 68‐4.

†The resistant parent Yuanza 9102.

‡The genomic regions determined based on the results of SP and RP.

§Chromosome.

¶SNP‐index of susceptible bulk – SNP index of resistant bulk.

††99% confidence interval upper side.

‡‡99% confidence interval lower side.

Therefore, two genomic regions were selected as the candidates for BWR, that is, 3.25–6.60 Mb on B02, and 13.55–16.40 Mb on B09 for further detailed analysis. Notably, the 3.35 Mb interval on B02 was significant (P < 0.01) irrespective of parent used as reference, while the 2.85 Mb interval on B09 was only significant (P < 0.01) when RP was used as reference.

Reconfirmation of identified genomic regions through genetic mapping

To validate and narrow down the genomic region identified for BWR on chromosome B02, KASP markers were successfully developed for 28 SNPs whose ΔSNP‐index were higher than statistical confidence of P < 0.01 (Table S3). Ten SNPs were located in the candidate genomic region while the remaining 18 SNPs were from nearby regions where less significant peaks of ΔSNP‐index were observed (Figure 3a–h). Genotyping of the RIL population was performed with the 28 KASP markers followed by genetic map construction with map length of 38.14 cM (Table S4). Composite interval mapping (CIM) confirmed the existence of the QTL and its stable in expression across three environments. This QTL was located between marker Araip_B02_3740746 and Araip_B02_5804063, with 14.4%–29.32% phenotypic variations explained (PVE) and LOD scores of 6.53–17.23 (Table 2). Single‐marker analysis (SMA) also confirmed the significant association between BWR and the three markers in the interval (Table S4). Therefore, the candidate genomic region on B02 was validated and narrowed down to a 2.07 Mb interval (3.74–5.81 Mb).

Table 2

Genetic and QTL mapping for bacterial wilt resistance using newly developed SNP markers on chromosome B02

Environment	Position (cM)	LOD value	Marker interval	Additive effect	PVE* (%)
Hongan2014	3.3	6.53	Araip_B02_3740746–Araip_B02_5804063	12.91	16.03
Nanchong2015	5.3	6.84	Araip_B02_3740746–Araip_B02_5804063	15.04	14.41
Nanchong2016	6.3	17.23	Araip_B02_3740746–Araip_B02_5804063	18.64	29.32

Phenotypic variation explained.

Similarly, eight KASP markers were successfully developed to validate the genomic region identified for BWR on chromosome B09 (Table S3). Genotyping of the RIL population was performed with these KASP markers, and a genetic linkage map was constructed with 4.86 cM length (Table S4). Unfortunately, the CIM analysis could not detect any QTL for BWR across three environments. In addition, SMA found that only marker Araip_B09_14990612 showed low significance in only Nanchong2016 environment with 5.42% PVE and 2.53 LOD score (Table S4). These results indicated that the 2.85 Mb interval on B09, which was identified only with the RP assembly as reference, did not harbour major and stable QTL for BWR.

In addition, the traditional QTL analysis was conducted using the same phenotyping data and the previously published SSR‐based genetic linkage map with 830 loci (Luo et al., 2017). CIM identified five additive QTLs with 4.63%–40.19% PVE across three environments (Table S5). Only one major QTL, qBWRB02.1, was identified on the chromosome B02 (13.02%–40.19% PVE), and it was stable in expression across three environments while no QTL was identified on B09. The flanking SSR markers of qBWRB02.1, that i Ad02A625 and Ai02B2356, were located at 0.81 and 5.95 Mb on the chromosome B02 of A. ipaensis (Figure S8). The above results confirmed that the 2.07 Mb genomic region (3.74–5.81 Mb) on B02 harboured a major and stable QTL for BWR and the resistant allele was from Yuanza 9102.

Putative candidate genes for bacterial wilt resistance

The genomic region spanning 2.07 Mb on chromosome B02 had 348 effective SNPs. Of the 348 SNPs, 246 SNPs were identified irrespective of parent used as reference; however, 76 SNPs were specifically identified when SP assembly was used as reference while 26 SNPs were specifically identified when RP assembly was used as reference (Tables S6‐S7). Function annotation analysis of the 348 SNPs found that 260 SNPs were intergenic, 20 intronic, 1 in 5′ UTR and 2 in 3′ UTR, 16 synonymous and 49 nonsynonymous (three resulted in stop codon). The 49 nonsynonymous SNPs affected 19 putative candidate genes for BWR (Table 3).

Table 3

Nonsynonymous SNPs in putative candidate genes in the genomic region for bacterial wilt resistance on chromosome B02

Gene	Position (bp)	SP base	RP base	SB base	RB base	Amino acid change	Function	Reference assembly
Araip.G1MIP	3740746	A	C	A	C	Phe789Leu	Disease resistance protein (CNL)	Both†
Araip.XX28B	3862695	T	A	T	A	Asp694Val	Extra‐large guanine nucleotide‐binding protein	SP‡
Araip.EK9PX	3889426	T	A	T	A	Glu15Val	UDP‐Glycosyltransferase superfamily protein	Both
Araip.VE4DY	3919494	T	G	T	G	Phe72Leu	Disease resistance protein (CNL)	Both
	3920191	C	G	C	G	His305Asp		Both
	3920337	A	C	A	C	Glu353Asp		SP
Araip.05JB8	3939153	C	G	C	G	Phe452Leu	Disease resistance protein (CNL)	SP
Araip.YCN52	4067048	C	T	C	T	Pro166Ser	Disease resistance protein (CNL)	SP
	4067057	T	G	T	G	Leu169Val		SP
	4067061	G	C	G	C	Cys170Ser		SP
	4067641	T	G	T	G	Met363Ile		RP§
	4067661	G	A	G	A	Lys370Arg		RP
	4068264	C	G	C	G	Thr571Arg		SP
	4068285	T	A	T	A	Phe578Tyr		SP
	4068291	T	A	T	A	Phe580Tyr		SP
	4068297	A	G	A	G	Glu582Gly		SP
	4068303	T	G	T	G	Phe584Cys		SP
	4068305	C	G	C	G	Pro585Ala		SP
	4068362	G	C	G	C	Gly604Arg		SP
	4068498	T	G	T	G	Phe649Cys		SP
	4068809	T	G	T	G	Tyr753Asp		SP
Araip.U0YME	4109028	G	C	G	C	Lys36Asn	Disease resistance protein (CNL)	SP
	4109068	G	A	G	A	Asp50Asn		RP
	4109096	G	A	G	A	Arg59His		RP
	4110068	G	T	G	T	Cys383Phe		Both
	4110086	C	G	C	G	Ser389*		SP
	4110521	A	G	A	G	Gln534Arg		Both
	4110629	T	C	T	C	Ile570Thr		Both
	4110719	A	T	A	T	Gln600Leu		SP
	4110739	G	A	G	A	Asp607Asn		SP
Araip.WF303	4123931	G	A	G	A	Val229Ile	Disease resistance protein (CNL)	RP
Araip.65IVT	4135837	T	G	T	G	Phe72Leu	Disease resistance protein (CNL)	SP
	4136087	C	G	C	G	Pro156Ala		SP
	4137246	C	T	C	T	Thr542Ile		RP
	4137389	A	C	A	C	Thr590Pro		SP
	4137449	G	C	G	C	Glu610Gln		SP
Araip.AA1CY	4548625	C	A	C	A	Glu556*	DYNAMIN‐like 1B	SP
Araip.ZD3PG	4726643	T	C	T	C	Ser133Leu	Auxin response factor	Both
Araip.R2VR0	4875274	C	G	C	G	Tyr759*	Vacuolar protein‐sorting protein bro1	Both
Araip.CZ54N	5447610	C	T	C	T	Pro168Ser	Serine/threonine‐protein phosphatase	Both
	5447943	A	G	A	G	Gly240Ser		Both
Araip.PKI4P	5635075	G	A	G	A	Gly31Arg	Pentatricopeptide repeat superfamily protein	Both
	5635294	T	C	T	C	Arg104Trp		Both
Araip.3CL97	5661423	G	C	G	C	His463Gln	Vacuolar protein‐sorting protein bro1	Both
Araip.HY224	5724170	G	C	G	C	Leu21Val	ALG‐2 interacting protein X‐like	Both
Araip.V2RJA	5778260	C	T	C	T	Asp447Asn	Vacuolar protein‐sorting protein bro1	Both
	5778267	A	C	A	C	Ile444Met		Both
Araip.E26FT	5787552	A	G	A	G	Gly5Arg	Unknown protein	Both
Araip.CMC8E	5804063	T	G	T	G	Ile193Met	Protein kinase family protein	Both

†The SNP was identified irrespective of parent used as reference.

‡The SNP was identified specifically with the susceptible parent Xuzhou 68‐4 as reference.

§The SNP was identified specifically with the resistant parent Yuanza 9102 as reference.

Notably, seven of the 19 putative candidate genes were predicted to code for the NBS‐LRR type disease resistance proteins, including Araip.G1MIP, Araip.VE4DY, Araip.05JB8, Araip.YCN52, Araip.U0YME, Araip.WF303 and Araip.65IVT (Figure 4). The putative candidate gene Araip.E26FT codes for protein with unknown function. The putative candidate gene Araip.PKI4P codes for pentatricopeptide repeat superfamily protein. Three putative candidate genes namely Araip.R2VR0, Araip.3CL97 and Araip.V2RJA code for vacuolar protein‐sorting protein bro1. The other putative candidate genes code for various kinds of proteins, including ALG‐2 interacting protein X‐like (Araip.HY224), auxin response factor (Araip.ZD3PG), calcium‐dependent protein kinase (Araip.CMC8E), DYNAMIN‐like 1B (Araip.AA1CY), extra‐large guanine nucleotide‐binding protein 1‐like (Araip.XX28B), serine/threonine‐protein phosphatase 7 long form homolog (Araip.CZ54N) and UDP‐Glycosyltransferase superfamily protein (Araip.EK9PX; Table 3).

Figure 4

The identified candidate NBS‐LRR resistance proteins associated with bacterial wilt resistance. CC: Coiled‐coil domain. NB: Nucleotide‐binding site domain. LRR: Leucine‐rich repeat domain. The positions of amino acid changes caused by nonsynonymous SNPs were shown in yellow colour.

Validation of diagnostic markers

To validate diagnostic markers for BWR, the three SNP markers located in the 2.07 Mb genomic region on B02 were validated initially on 21 resistant and nine susceptible breeding lines derived from the crossing between Yuanza 9102 and a susceptible cultivar Zhonghua 5. Of the three markers, two markers (Araip_B02_3740746 and Araip_B02_5804063) showed clear differentiation between resistant and susceptible breeding lines. The first diagnostic marker, Araip_B02_3740746, amplified the CC alleles in the resistant parent while the AA alleles in the susceptible parent. The second diagnostic marker, Araip_B02_5804063, amplified the GG alleles in the resistant parent while the TT alleles in the susceptible parent (Figure 5a). All 21 resistant breeding lines showed CCGG alleles while the nine susceptible lines showed AATT alleles (Figure 5b, Table S8).

Figure 5

The validated diagnostic markers for bacterial wilt resistance. (a) The primers of the two diagnostic markers. †Sequences of the allele‐specific primers do not include the tail sequences that interact with the fluor‐labelled oligos in the KASP reaction mix. (b) Scatter plot of the survival rates (%) of breeding lines with different genotypes using the diagnostic markers. (c) Box plots for comparing the difference of survival rates (%) between RILs with different genotypes. Box plots were generated with the ggpubr package in R software. In each box, centre line shows the median; box limits indicate the 25th and 75th percentiles; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles. **** indicated the difference is significant at 0.0001 level.

Secondly, two resistant cultivars, Yuanza 5 and Zhonghua 21, were derived from the crossings Yuanza 9102 × Shiyouhong 4 and Yuanza 9102 × Zhonghua 5, respectively. The diagnostic markers amplified the resistant alleles CCGG in both of them; thus, their BWR might come from Yuanza 9102 (Table S8). In addition, the diagnostic markers were genotyped in a panel consisting of 13 susceptible and 13 other resistant cultivars (Table S8). All of them amplified the AATT alleles. Therefore, the resistant allele from Yuanza 9102 might be different from the 13 resistant cultivars.

Thirdly, RILs of the mapping population were classified into two classes, CCGG and AATT, according to the genotypes of the two diagnostic markers. The mean survival rate of the 59 CCGG RILs was 74.45% and significantly higher than that of the 78 AATT RILs (34.17%), indicating that the resistant allele from Yuanza 9102 significantly improved BWR in all three environments (~40% survival rates on average, Figure 5c, Table S9).

The above results indicated that markers, Araip_B02_3740746 and Araip_B02_5804063, would be useful for selecting breeding lines with resistance to bacterial wilt.

Discussion

The recently completed genome sequences of the diploid ancestors of cultivated peanut (Bertioli et al., 2016) enabled the utilization of NGS‐based technologies in rapid trait mapping in peanut (Pandey et al., 2016). The QTL‐seq approach only requires WGRS data of extreme bulks and parents of mapping population, making it cost‐effective and fast (Takagi et al., 2013). This approach has been proved to be very successful with RIL populations in several crops including rice (Takagi et al., 2013), peanut (Pandey et al., 2017), cucumber (Wei et al., 2016), chickpea (Deokar et al., 2018) and pigeonpea (Singh et al., 2016). The present study successfully applied the QTL‐seq approach (Figure 2) to identify genomic regions and candidate genes for BWR using the sequencing data of parental genotypes and pooled samples of RIL population (Yuanza 9102 × Xuzhou 68‐4).

In order to increase the reliability of identified genomic regions and candidate genes, we used resistant as well as susceptible parent separately as reference parent in the QTL‐seq approach. According to the original QTL‐seq approach (Takagi et al., 2013), only one parent was required to be used as reference to analyse the reads of two extreme bulks, and previous reports followed this principle in multiple crops (Deokar et al., 2018; Hua et al., 2016; Pandey et al., 2017; Singh et al., 2016). Recently, in our one such study for shelling percentage in peanut, we found that when both parents were sequenced and used to analyse the reads of two extreme bulks in parallel, the comparative analysis of the two sets of results generated by QTL‐seq increased the mapping efficiency and accuracy of shelling percentage in peanut (Luo et al., 2018). Therefore, we used this similar strategy in the present study for mapping BWR (Figure 2) and found significance differences between the two sets of results as well. Interestingly, the genomic region on B09 was identified only when RP assembly was used as reference and was proved to be a minor and unstable QTL by SMA. By contrast, the genomic region on B02 was identified irrespective of parents and was validated to be a major and stable QTL by both CIM and SMA. However, the positions of candidate SNPs associated with BWR in the 2.07 Mb interval on B02 were not the same, that is, 76 SNPs were specifically identified when SP assembly used as reference while 26 SNPs were specifically identified when RP assembly used as reference. The above results indicated that the choice of different parent in the QTL‐seq approach would affect the results. The large differences of the two sets of results might due to different level of alignment errors upon using different parent‐based reference genomes. For example, the numbers of reads mapped to the locations of identified SNPs varied a lot when different parents were used for sequence analysis of extreme bulks (Figure S9). The large alignment differences of reads might due to the complexity of polyploid genomes (e.g. the similarity between A and B subgenomes). Irrespective of the reasons for the large difference, the pipeline used in the present study (Figure 2), through comparing and merging the two sets of results generated by both parents, provided precise and accurate discovery of candidate genomic regions, genes and SNPs markers for BWR.

The QTL‐seq approach could be successful in identifying candidate regions for major and stable QTLs of BWR but not for minor QTLs. Through the classical QTL mapping with previously reported SSR marker‐based genetic linkage map (Luo et al., 2017), one major QTL qBWRB02.1 and four minor QTLs were identified for BWR. The QTL qBWRB02.1 might span a 5.14 Mb (0.81–5.95 Mb) interval on chromosome B02 based on its flanking SSR markers. Using the QTL‐seq approach (Figure 2), a 2.07 Mb (3.74–5.81 Mb) genomic region on B02 was identified for BWR across three environments (Table 2) and found to be co‐localized with the major QTL qBWRB02.1 (Figure S8). However, no corresponding genomic regions were identified for the other four minor QTLs, owing to the fact that their effects were relative low (4.63%–6.48% PVE) and their expressions were not stable across environments. The extreme bulks of QTL‐seq were constructed based phenotyping data of multiple environments that might be one of the reasons for missing to detect minor and unstable QTLs. Wang et al. (2018b) identified the same major QTL as well as three minor QTLs on B02 based on the phenotype of 188 RILs and genetic map constructed by double‐digest restriction‐site‐associated DNA sequencing (ddRADseq). However, the density of SNPs identified by ddRADseq was low, which might be limited by the double‐digest restriction‐sites, and provided limited information for the identification of genomic regions and candidate genes (Wang et al., 2018b). In QTL‐seq, high dense SNPs were identified by resequencing in the genomic region of qBWRB02.1, providing key information for discovering candidate genes and diagnostic markers. These results illustrated that the genomic region identified for the major QTL qBWRB02.1 with the QTL‐seq pipeline (Figure 2) was robust and precise.

The qBWRB02.1 on chromosome B02 would be a novel and valuable resistant allele in breeding elite cultivars with enhanced BWR. Zhao et al. (2016) identified a major QTL qBW‐1, whose resistance was from peanut cultivar Yueyou 92. The genomic location of qBW‐1 was not reported; however, the linked marker, SNP79, was reported to be located in a BAC clone (GenBank Accession Number HQ637177.1). Through BLAST, we found this BAC clone located on chromosome B04 (Table S10). In addition, Zhang et al. (2017) isolated NBS‐LRR gene AhRRS5 in Yueyou 92 and reported the overexpression of AhRRS5 in tobacco enhanced BWR. Through BLAST, we found AhRRS5 located on chromosome A05 or B05 (Table S10). Based on the pedigree tracking, the source of resistance in Yueyou 92 is from a Chinese landrace Xiekangqing while the source of resistance in Yuanza 9102 is from wild species (Shen et al., 2015). Therefore, we assume that the qBWRB02.1 from Yuanza 9102 would be different from the reported loci or genes from Yueyou 92. Moreover, we found that the resistant allele of qBWRB02.1 from Yuanza 9102 was successfully transferred during the breeding of resistant cultivars Yuanza 5 and Zhonghua 21 but was not present in the other cultivars (Table S8).

Wild Arachis species might be important source to enhance the diversity of resistant genes in peanut breeding. The peanut cultivar Yuanza 9102 was developed through interspecific hybridization between the cultivated peanut Baisha 1016 and A. diogoi (Shen et al., 2015). Since Baisha 1016 is susceptible to bacterial wilt, the resistance of Yuanza 9102 might be from A. diogoi. Although the resistant allele is from A. diogoi, it has no negative genetic effect on yield trait and could be useful in breeding of high BWR and yield cultivars, e.g. Zhonghua 21 (4439 kg/ha). The genus Arachis contains 79 wild species and provides vast genetic diversity (Krapovickas et al., 2007). These wild species harbour new and diverse sources of R‐genes which could be transferred into the cultivated gene pool by multiple pre‐breeding strategies (Moore, 2015). For example, the wild Arachis species have been harnessed to improve the resistance to rust and late leaf spot in peanut (Sharma et al., 2017). Therefore, the pyramiding of the resistance QTLs or genes from both cultivated peanut and its wild relatives through GAB might provide higher level or more stable resistance to the bacterial wilt disease.

The QTL‐seq pipeline was effective in the identification of candidate SNPs associated with BWR, which could be developed into genetic markers based on different genotyping techniques for fine mapping or using as diagnostic markers after validation. In the present study, 348 SNPs were identified in the 2.07 Mb genomic region on chromosome B02, accounting ~1‰ of the total number (332 822) of genome‐wide SNPs. Of the 348 SNPs, 25.29% were located in genic region, providing information on putative candidate genes for BWR (Hartwig et al., 2012). There were 49 nonsynonymous SNPs affected 19 putative candidate genes for BWR (Table 3). In the present study, SNPs were successfully transformed into KASP markers, because genotyping with KASP markers is cheaper than SSR markers (Steele et al., 2018). We found that the KASP markers, Araip_B02_3740746 and Araip_B02_5804063, showed clear differentiation between resistant and susceptible genotypes in a diverse validation panel consisting of breeding lines and cultivars (Table S8). Using the two markers as diagnostic markers, 59 RILs were detected to possess resistant alleles from Yuanza 9102, and they improved ~40% survival rates on average when compared to the 78 RILs with susceptible alleles from Xuzhou 68‐4 (Table S9). These results indicated the potential of these KASP markers in tracking the resistant alleles for BWR in the breeding programmes.

R‐genes play a key role in the plant immune system (Jones and Dangl, 2006; Spoel and Dong, 2012). The proteins encoded by R‐genes usually share common domains such as nucleotide‐binding site (NBS), Toll interleukin receptor‐like (TIR), leucine‐rich repeat (LRR) and kinase domain (Sanseverino et al., 2010). Most of R‐genes cloned encode NBS‐LRR proteins, which could be divided into TNL and CNL based on the presence or absence of TIR domains. Receptor‐like proteins (RLP) and the receptor‐like kinases (RLK) are two major transmembrane receptor proteins encoded by R‐genes. In A. thaliana, the RLK (Godiard et al., 2003) and TNL (Deslandes et al., 2002) genes have been proven to control BWR. Of the 19 putative candidate genes identified in the present study, seven were CNL type R‐genes, including Araip.G1MIP, Araip.VE4DY, Araip.05JB8, Araip.YCN52, Araip.U0YME, Araip.WF303 and Araip.65IVT (Figure 4). Notably, maximum number of nonsynonymous SNPs (14) was identified in the Araip.YCN52 gene, and one nonsynonymous SNPs in the Araip.U0YME gene resulted in stop codon. Recently, Zhang et al. (2017) overexpressed the peanut CNL gene, AhRRS5, in tobacco and found that the BWR of transgenic lines enhanced. Therefore, the CNL genes might be an important type of R‐genes for BWR in peanut.

The remaining 12 putative candidate genes might be involved in the BWR as well. The gene Araip.PKI4P might encode for pentatricopeptide repeat superfamily protein whose dysfunction rendered A. thaliana more susceptible to pathogenic bacteria (Park et al., 2014). The gene Araip.ZD3PG was predicted to code for auxin response factor, which may involve in the regulation of auxin pathway. Recent reports showed that the repression of auxin pathway improved the resistance to bacterial pathogen (French et al., 2018; Navarro, 2016). The gene Araip.CMC8E might code for calcium‐dependent protein kinase, which could perceive changes in intracellular calcium concentrations to transmit defence signal in response to PAMPs and pathogen effectors (Liu et al., 2017; Wang et al., 2018a). The gene Araip.CZ54N was predicted to encode for serine/threonine‐protein phosphatase, which was reported to be involved in hypersensitive response (Zhou et al., 1995). The two G proteins encoded by Araip.AA1CY and Araip.XX28B might be involved in the signal transduction that resulted in a change in cell function (Zhu et al., 2009). The gene Araip.EK9PX was predicted to encode UDP‐Glycosyltransferase superfamily protein, whose reduced expression improved resistance to bacteria in A. thaliana (Park et al., 2011). Based on these findings, these genes together with the seven R‐genes should be targeted as candidates for fine mapping and function validation.

In conclusion, the QTL‐seq approach is a powerful method and successfully identified genomic regions and candidate genes for major and stable QTLs for BWR. The present study identified a 2.07 Mb genomic region on chromosome B02 for BWR containing seven R‐genes. Most importantly, the diagnostic markers developed and validated in this study could be used in GAB in future breeding programmes to accelerate development of peanut varieties with enhanced BWR.

Materials and methods

Plant materials

The peanut cultivar Yuanza 9102, developed by Henan Academy of Agricultural Sciences, China in 2002 through interspecific hybridization between a cultivated peanut cultivar Baisha 1016 and a wild specie A. diogoi, was stably resistant to bacterial wilt during multiple years of field assessment (Figure 1a and b). Yuanza 9102 was crossed with a susceptible cultivar Xuzhou 68‐4 and advanced to the F9 generation following single seed decent (SSD) method (Luo et al., 2018). A total of 195 RILs were used for trait mapping of BWR. To validate the diagnostic markers, a panel consisting of 21 resistant and nine susceptible breeding lines, which derived from the crossing between Yuanza 9102 and another susceptible cultivar Zhonghua 5, and another panel consisting of 13 susceptible cultivars and 13 resistant cultivars were used for genotyping with newly developed allele‐specific markers. In addition, two resistant cultivars, Yuanza 5 and Zhonghua 21, were included in this study. Yuanza 5 derived from the cross between Yuanza 9102 and a susceptible cultivar Shiyouhong 4 by Henan Academy of Agricultural Sciences in 2013. Zhonghua 21 derived from the cross between Yuanza 9102 and a susceptible cultivar Zhonghua 5 by Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences in 2012.

Phenotyping of bacterial wilt resistance

In order to assess the resistance of peanut materials to bacterial wilt, a disease nursery was constructed in the Hongan city, Hubei province, China (31.36° N, 114.61° E) in 1964, and another disease nursery was constructed in the Nanchong city, Sichuan province, China (30.67° N, 106.06° E) in 2007. The BWR of peanut materials were screened in a randomized block design with three replications in each disease nursery following standard procedures (Jiang et al., 2006; Wang et al., 2018b). The plants killed by the bacterial wilt disease were recorded until harvesting, and the survival rate was calculated with the formula:

The RIL population (Yuanza 9102 × Xuzhou 68‐4) was planted in the disease nursery in Hongan during 2014 and 2015 as well as in Nanchong during 2015 and 2016, which were designated as four environments: Hongan2014, Hongan2015, Nanchong2015 and Nanchong2016. The Hongan2015 environment was excluded in the present study due to its inadequate disease incidence. The survival rates of 30 breeding lines used to validate diagnostic markers were successfully screened in Hongan during 2004 and 2006. The survival rates of 26 cultivars used to validate diagnostic markers were assessed in Hongan during 2011, 2012 and 2013.

Construction of extreme bulks and Illumina sequencing

Unexpanded leaflets of the RIL population in F9 generation were collected and used to extract high‐quality genomic DNA with a modified CTAB method (Doyle, 1990). The mean survival rate of each RIL was calculated based on the precise phenotyping data obtained in Hongan2014, Nanchong2015 and Nanchong2016. To develop the resistant bulk (RB) for BWR, the same amounts of DNA from 15 RILs with high mean survival rates were pooled. Similarly, the same amounts of DNA from 15 RILs with low mean survival rates were pooled to construct the susceptible bulk (SB) for BWR.

The genomic DNA of the extreme bulks (RB and SB) and the parents (Yuanza 9102 and Xuzhou 68‐4) was used to construct Illumina libraries (~350 bp) following the protocol of the NEBNext Ultra II DNA Library Prep Kit for Illumina. Pair‐ended reads with 150 bp length of the four libraries were generated using the Illumina NovaSeq platform with NovaSeq 6000 S4 Reagent Kit. The raw sequencing data have been submitted in the NCBI Sequence Read Archive (SRA) under BioProject accession number PRJNA528382.

Identification of genomic regions for bacterial wilt resistance

To map the QTLs for BWR, a pipeline based on the QTL‐seq approach (Takagi et al., 2013), as shown in Figure 2, was used in the present study. The high‐quality reads with more than 90% of nucleotides showed Phred quality scores higher than 30 were extracted from the raw sequencing reads and used for following analysis.

To develop reference‐guided assembly of parents, the genome assemblies of A. duranensis and A. ipaensis were downloaded from the PeanutBase website (http://www.peanutbase.org; Bertioli et al., 2016) and used as the reference. For developing a reference‐guided assembly for the susceptible parent Xuzhou 68‐4 (the SP assembly), the high‐quality reads of SP were aligned to the reference using the BWA software (Li and Durbin, 2009) and further refined with the Coval software (Kosugi et al., 2013). Then, SNPs between the SP and reference were called with the Samtools software (Li et al., 2009) and filtered with the Coval scripts (Kosugi et al., 2013). Finally, the SP assembly was developed by substituting the reference bases with alternative bases at the positions of filtered SNPs. Similarly, a reference‐guided assembly for the resistant parent Yuanza 9102 (the RP assembly) was constructed.

To identify the genomic regions controlling BWR, the high‐quality reads from SB and RB were further filtered by removing unpaired reads. To equalize the number of reads between bulks, the filtered reads of RB were randomly reduced to the same amount of the filtered reads of SB. Using the SP assembly as reference, the equalized reads of SB and RB were aligned and SNPs were called with default parameters of the QTL‐seq_framework1.4.4 developed by the Iwate Biotechnology Research Center, Japan (Takagi et al., 2013). For each SNP in each bulk, SNP‐index was calculated with the formula: (Figure S1). By subtracting SNP‐index of RB from SNP‐index of SB, ΔSNP‐index for each SNP was calculated. The SNPs with SNP‐index <0.3 or read depth <10 in both of the bulks were filtered out in followed analysis. Then, sliding window analyses with 2 Mb interval and 50 kb increment were conducted for SNP‐index and ΔSNP‐index of the filtered SNPs ([Link], [Link], [Link], [Link], [Link], [Link]). The genomic regions, whose ΔSNP‐index significantly deviated from the statistical confidence intervals under the null hypothesis of no QTLs at P < 0.01 level and SNP‐index significantly deviated from 0.5 in both bulks, were selected as candidates controlling BWR. On the other hand, the equalized reads of SB and RB were aligned to the RP assembly, and the second set of genomic regions for BWR was identified with the same approach. Finally, the two set of genomic regions identified with the SP and RP assemblies, respectively, were combined and considered as the genomic regions for BWR.

SNP marker development and validation of candidate genomic regions identified

To validate the genomic regions for BWR identified above, SNPs with SNP‐index significantly deviated from 0.5 in both bulks and ΔSNP‐index significantly higher or lower than the statistical confidence intervals under the null hypothesis of no QTLs at P < 0.01 level were selected for marker development. The 100 bp upstream and downstream sequences of the selected SNPs were extracted from the genomic sequences of A. ipaensis and blast to the genomic sequences of A. duranensis and A. ipaensis. Only SNPs with specific hit were used as template to developed KASP markers (Semagn et al., 2014). Each KASP marker contained two allele‐specific forward primers and one common reverse primer. SNP markers were validated with the parents of the RIL population (Table S4). RILs of the mapping population were genotyped with the validated SNP markers. The JoinMap 4.0 software (Van Ooijen, 2006) was used to construct the genetic map for targeted linkage group/chromosome. Genetic distances were estimated with the Kosambi mapping function (Kosambi, 1943). The survival rates of the RIL population generated in three environments, that is Hongan2014, Nanchong2015 and Nanchong2016, were used for conducting the QTL analysis. The QTLs for BWR were identified by both composite interval mapping and single‐marker analysis using the Windows QTL Cartographer 2.5 software (Wang et al., 2012).

Identification of putative candidate genes

Of the SNPs identified using the SP assembly, effective SNPs that might be associated with BWR were identified with criteria: read depth of ≥10, SNP‐index lower than 0.5 in SB but higher than 0.5 in RB, and ΔSNP‐index significantly lower than the statistical confidence intervals lower side under the null hypothesis of no QTLs at P < 0.01 level (Table S6). On the other hand, effective SNPs were also screened from the SNPs identified using the RP assembly with criteria: read depth of ≥10, SNP‐index significantly higher than 0.5 in SB but lower than 0.5 in RB, and ΔSNP‐index significantly higher than the statistical confidence intervals upper side under the null hypothesis of no QTLs at P < 0.01 level (Table S7). Then, these effective SNPs were merged and function annotation was conducted with the SnpEff v3.0 software (Cingolani et al., 2012). The putative genes affected by nonsynonymous SNPs were considered as the candidate genes responsible for the BWR. The protein domains encoded by the putative candidate genes were predicted with the InterPro online server (Finn et al., 2017).

Conflict of Interest

The authors declare that there is no conflict of interest.

Author contributions

HL, MKP, BL, RKV and HJ conceived, designed and supervised the experiments. XR, YL, BL and HJ developed the RIL population. HL, XR, BW, JG, LH, XZ, YC, WC and NL conducted field trials and phenotyping. HL, BW, JG and WC performed DNA extraction and genotyping. HL, MKP, AWK, BL, RKV and HJ performed the QTL‐seq analysis and interpreted the results. HL prepared the first draft and HL, MKP, YL, BL, RKV and HJ contributed to the final editing of manuscript. All authors read and approved the final manuscript.

Figure S1 Alignment, SNP identification and calculation of SNP‐index and ΔSNP‐index for bacterial wilt resistance.

Click here for additional data file.

Figure S2 Genome‐wide SNP‐index plots of susceptible bulk with the susceptible parent Xuzhou 68‐4 as reference.

Click here for additional data file.

Figure S3 Genome‐wide SNP‐index plots of resistant bulk with the susceptible parent Xuzhou 68‐4 as reference.

Click here for additional data file.

Figure S4 The ΔSNP‐index plot obtained by subtraction of resistant bulk SNP‐index from susceptible bulk SNP‐index using the susceptible parent Xuzhou 68‐4 as reference.

Click here for additional data file.

Figure S5 Genome‐wide SNP‐index plots of susceptible bulk with the resistant parent Yuanza 9102 as reference.

Click here for additional data file.

Figure S6 Genome‐wide SNP‐index plots of resistant bulk with the resistant parent Yuanza 9102 as reference.

Click here for additional data file.

Figure S7 The ΔSNP‐index plot obtained by subtraction of resistant bulk SNP‐index from susceptible bulk SNP‐index using the resistant parent Yuanza 9102 as reference.

Click here for additional data file.

Figure S8 Co‐localization of QTLs identified by SNP‐ and SSR‐based genetic map for bacterial wilt resistance.

Click here for additional data file.

Figure S9 The scatterplots of depth of reads mapped to the identified SNPs in extreme bulks using different parent as reference.

Click here for additional data file.

Table S1 Details on whole genome re‐sequencing data.

Table S2 Chromosome‐wise SNPs distribution between extreme bulks for bacterial wilt resistance.

Table S3 Details of newly developed SNP markers.

Table S4 Single‐marker analysis for bacterial wilt resistance using validated SNP markers.

Table S5 QTLs identified for bacterial wilt resistance using the SSR markers‐based genetic map.

Table S6 Effective SNPs identified for bacterial wilt resistance using the susceptible parent Xuzhou 68‐4 assembly as reference.

Table S7 Effective SNPs identified for bacterial wilt resistance using the resistant parent Yuanza 9102 assembly as reference.

Table S8 The genotypes of diagnostic markers in the breeding lines and cultivars with different resistance to bacterial wilt.

Table S9 Phenotypic effect of the major and stable QTL on chromosome B02 for bacterial wilt resistance in the RIL population.

Table S10 Blast results of the BAC clone and gene associated with bacterial wilt resistance in Yueyou 92.

Click here for additional data file.