A Major and Stable QTL for Bacterial Wilt Resistance on Chromosome B02 Identified Using a High-Density SNP-Based Genetic Linkage Map in Cultivated Peanut Yuanza 9102 Derived Population

Bacterial wilt (BW) is one of the important diseases limiting the production of peanut (Arachis hypogaea L.) worldwide. The sufficient precise information on the quantitative trait loci (QTL) for BW resistance is essential for facilitating gene mining and applying in molecular breeding. Cultivar Yuanza 9102 is BW resistant, bred from wide cross between cultivated peanut Baisha 1016 and a wild diploid peanut species A. chacoense with BW resistance. In this study, we aim to map the major QTLs related to BW-resistance in Yuanza 9102. A high density SNP-based genetic linkage map was constructed through double-digest restriction-site-associated DNA sequencing (ddRADseq) technique based on Yuanza 9102 derived recombinant inbred lines (RILs) population. The map contained 2,187 SNP markers distributed on 20 linkage groups (LGs) spanning 1566.10 cM, and showed good synteny with AA genome from A. duranensis and BB genome from A. ipaensis. Phenotypic frequencies of BW resistance among RIL population showed two-peak distribution in four environments. Four QTLs explaining 5.49 to 23.22% phenotypic variance were identified to be all located on chromosome B02. The major QTL, qBWB02.1 (12.17–23.33% phenotypic variation explained), was detected in three environments showing consistent and stable expression. Furthermore, there was positive additive effect among these major and minor QTLs. The major QTL region was mapped to a region covering 2.3 Mb of the pseudomolecule B02 of A. ipaensis which resides in 21 nucleotide-binding site -leucine-rich repeat (NBS-LRR) encoding genes. The result of the major stable QTL (qBWB02.1) not only offers good foundation for discovery of BW resistant gene but also provide opportunity for deployment of the QTL in marker-assisted breeding in peanut.


INTRODUCTION
Cultivated peanut (Arachis hypogaea L.) is an important oilseed crop in the world (Bertioli et al., 2016;Varshney et al., 2017). The bacterial wilt (BW) caused by Ralstonia solanacearum is the most important soil-borne bacterial disease in peanut, especially destructive in China. It usually causes 10-30% yield loss. In the extreme cases, it results in yield loss more than 50% even extinction (Wicker et al., 2007;Yu et al., 2011). This disease is also seriously harmful to tomato, eggplant, chili, potato, tobacco, and other crops (Cao et al., 2009;Liao et al., 2010;Prasath et al., 2011;Narusaka et al., 2014;Xiao et al., 2015). However, the efficacy of chemical control, agricultural cultivation of crop rotation and biological control of BW is very limited, due to the environmental pollution resulted from chemical prevention or the higher cost for rotation and biological control. The most effective and economic measure is to breed and plant BW resistant cultivar (Lyu et al., 2010;Sunkara et al., 2014;Bhatnagar-Mathur et al., 2015).
The quantitative trait loci (QTL) controlling target trait is the important basis for marker-assisted selection (MAS) (Varshney et al., 2005(Varshney et al., , 2006. Many efforts have been made to identify molecular markers linked to the BW resistance for molecular breeding in peanut (Jiang et al., 2007;Ren et al., 2008;Peng et al., 2010;Hong et al., 2011;Zhao et al., 2016). Ren et al. (2008) identified two markers linked to BW resistance with genetic distance of 8.12 and 11.46 cM using 45 polymorphic amplified fragment length polymorphism (AFLP) markers and bulk segregant analysis (BSA) analysis. Peng et al. (2010) detected three QTLs for BW resistance based on a linkage map consisting of 98 AFLP markers and covering 285 cM map length. The analysis of these studies was based on the linkage maps with insufficient markers that linked with the resistant trait, thus limiting the use of MAS in the BW resistance breeding. Therefore, it is necessary to construct dense genetic map for precise QTL location.
Furthermore, identifying and pyramiding of QTLs derived from different resistant source is helpful to improve the resistant level and stability in BW resistance cultivar (Shi et al., 2008). Recently, two major QTLs based on a Yueyou 92 derived population was reported (Zhao et al., 2016). The source of resistance in Yueyou 92 can be traced back to Xiekangqing, a cultivated peanut with high level BW resistance. Meanwhile, Yuanza 9102, another popular cultivated peanut, has obvious resistance to BW in different condition too. However, the BW resistance of Yuanza 9102 is contributed by its parental line A. chacoense, a diploid wild peanut which has rich resistant resource. Whether the genetic basis of BW resistance in Yuanza 9102 and Yueyou 92 are same is unclear. It is worthy of understanding the genetic basis of BW resistance in Yuanza 9102 through QTL mapping.
With the development of the next generation sequencing (NGS), the high-throughput technology provides a convenient and quick platform for calling SNP markers and genotyping on a large scale (Baird et al., 2008;Peterson et al., 2012).
The double-digest restriction-site-associated DNA sequencing (ddRADseq) based on the next-generation sequencing has many advantages, such as reduced genomic complexity, high effectiveness, and moderate-cost (Baird et al., 2008;Davey et al., 2011Davey et al., , 2013Peterson et al., 2012;Gupta et al., 2015). Therefore, it has a wide range of applications on high-resolution genetic map, genome sequencing assembly as well as genetic analysis of complex genomic research (Peterson et al., 2012;Chen et al., 2013;Recknagel et al., 2013).
Molecular MAS provides an efficient tool for breeding as it offers rapid and precise selection of the targeted trait. In peanut, the technology has been used to improve several traits, such as rust resistance, nematode resistance, and the high oleic acid (Varshney et al., 2010(Varshney et al., , 2014aChu et al., 2011;Sukruth et al., 2015;Janila et al., 2016). However, the previous method for identification BW resistance mainly relies on phenotypic investigation in peanut field. Since the degree of the BW disease can be affected by many factors, such as soil conditions, bacteria concentration, rainfall and temperature, phenotypic investigation needs to be repeated for many different conditions. It is necessary to establish a quick, simple, and convenient method to improve the efficiency of peanut BW resistance selection. Combing precise QTL and closely linked molecular markers, resistant materials can be identified in the early generations using MAS technique. This will greatly improve the efficiency of peanut BW resistance breeding.
In this study, we employed ddRAD-seq approach to identify large scale genome-wide SNPs using the recombinant inbred line (RIL) population of Yuanza 9102 × Xuzhou 68-4. A SNPbased genetic linkage map was constructed and the syntenic analysis with AA and BB genome of diploid wild peanut was performed. Moreover, we identified QTLs for BW resistance in peanut and predicted candidate BW resistance genes for further investigation.

Plant Materials
A RIL population including 188 F 7 lines derived from the cross between Yuanza 9102 and Xuzhou 68-4 was used in this study. Yuanza 9102 is a popular cultivar in China that is resistant to R. solanacearum, while Xuzhou 68-4 is susceptible to R. solanacearum. Yuanza 9102 had wild-derived source of BW resistance which was derived from the cross of Baisha 1016 × A. chacoense (Xiong et al., 2003). The RIL lines were grown using randomized complete block design with three replications at two experimental locations. One was at experimental station of Nanchong (NC), Sichuan, China and the other was at experimental station of Hongan (HA), Hubei, China. Each accession was planted in three rows. Resistant parent Yuanza 9102 and susceptible parent Xuzhou 68-4 were planted after every 50 accessions as resistant and susceptible controls. There were 10 plants in each row with plant-to-plant distance of 10 cm. The row-to-row space was 30 cm apart. Genomic DNA was extracted from young leaf tissue essentially as described by Grattapaglia and Sederoff (1994).

ddRAD Library Construction and Sequencing
Genomic DNA was digested with EcoRI and MseI at 37 • C for 6 h and then heat inactivated at 65 • C for 90 min. The ligation with barcode adapter and common adapter, purification of pooled DNA template with specific adapters, and PCR enrichment were performed using a protocol described by Chen et al. (2013). Fragments of 250-450 bp were collected and dissolved in elution buffer for sequencing. The ddRAD library was sequenced on a Hiseq4000 next-generation sequencing platform. Raw Illumina reads with poor quality, contaminant sequences were filtered out using NGS QC Toolkit and data of individuals were separated according to barcode. The sequences of barcode were then removed and the clean data of each accession were used for subsequent analysis.

SNP Discovery and Genotyping
Using stacks-0.9998 software (Catchen et al., , 2013, the reads of the specific locus EcoRI and MseI flanking genomic DNA fragment were proceeded to cluster analysis and SNPs detection. The parameter of parents was set as: stacks-0.9998/ustacks -t gzfastq -r -d -m 5-M 2 -p 15, and the parameter of offspring was set as: stacks-0.9998/ustacks -t gzfastq -r -d -m 3 -M 2 -p 15. Using cstack software, all the loci of two parents were clustered as a catalog. The parameters was set as: stacks-0.9998/cstacks -b 1 -s QT0859 -s QT0860 -p 15 -n 3. Using sstack software, the stacks of RIL individuals were aligned to the catalog of parents and the parental alleles were determined in the progeny. The parameter was set as: stacks-0.99998/sstacksp 15 -b 1 -c * /result/batch_1 -s * /result/sample.fq -o * /result. Only the homozygous and polymorphic SNP markers were chosen for following analysis. In order to use the data more effectively, probabilistic PCA (PP) method was chosen to supplement the genotype data of each line (Fu, 2014). This method is based on PCA algorithm and has rapid computing speed and higher accuracy. Then, the proportion of missing data should be no more than 30%. The sequences of obtained SNPs were aligned against AA and BB genome of wild peanut A. duranensis and A. ipaensis (Bertioli et al., 2016), and only the sequences that uniquely mapped the genomes with the similarity and coverage over 95% were chosen. The remained SNPs were used for map construction. The developed SNP markers were designated as AHES (Arachis hypogea EcoRIdigestion genomic SNP) and AHMS (Arachis hypogea MseIdigestion genomic SNP) as described in Shirasawa et al. (2012).

Genetic Map Construction and Synteny Analysis
The genotypes of the 188 RIL F 7 individuals were used to construct a genetic map. LG assignments and marker order were based on the regression mapping of Joinmap 4.0 software (Van Ooijen, 2006). The map length was estimated using Haldane's mapping function. The graphic of the genetic map was represented using Mapchart 2.0 software (Voorrips, 2002). Markers were tested for segregation distortion by the chi-square test. The tags including the SNP markers of the genetic map were aligned with AA and BB genomes using blat software. The syntenic relationship was determined according to the physical location and genetic map distance and the graphic was displayed with circos figure.

Evaluation of Resistance to Bacterial Wilt
All individuals of the RIL population and their parents were planted in the BW disease nursery in Hongan and Nanchong. These disease nurseries were specially used for the evaluation of resistance to BW caused by R. solanacearum since 2007. The phenotype of BW resistance was investigated at Hongan in 2014 and 2015 and at Nanchong in 2015 and 2016. The method of evaluation of BW resistance was according to that described previously by Jiang et al. (2006) with small modification. At the seedling stage, flowering-pegging stage and harvest stage, the survival rate of each accession were investigated. Survival rate of each accession equals healthy plants/total plants × 100%. The resistance to BW was evaluated using the following standard: resistant, survival rate was of >80%; moderately resistant, survival rate ranged from 65 to 80%; moderate susceptible, survival rate ranged from 35 to 65%; susceptible, survival rate was <35%.

Statistical Analysis of Phenotyping Data
The broad-sense heritability was calculated as h 2 = σ 2 g /(σ 2 g + σ 2 ge / n + σ 2 e /nr) according to Hallauer and Miranda's method (Hallauer and Miranda, 1998). In the formula, σ 2 g represents the genetic variance, σ 2 ge represents the interaction variance of the genotype with environment, σ 2 e represents the variance of residual error, n represents the number of environments and r represents the number of replications. The calculation of the variance components was obtained using the SAS software by general linear model (GLM) procedure.

QTL Mapping for BW Resistance
Quantitative trait loci analyses were performed using the software WinQTLCart 2.5 (Wang S. et al., 2012) and program of composite interval mapping (CIM). The CIM analysis was performed using Model 6 and forward regression method. The number of control markers, window size, and walk speed were set as 5, 10, and 2 cM, respectively. The threshold of LOD for declaring the presence of a QTL was determined by 1000 permutation tests at a significance level of P < 0.01. QTLs had phenotypic variation explained more than 10% were considered as major QTLs. QTLs detected in different year/season were designated as consistent QTLs. The same QTL appeared in different locations were refereed as stable QTLs (Varshney et al., 2014b). QTLs were designated with an initial letter 'q' followed by the trait name and the LG corresponding chromosome, according to the previous described nomenclature (McCouch et al., 1997). If more than one QTL was detected for the same trait and linkage group (LG), a number was added after the LG. A graphic representation of the QTL was generated using Mapchart 2.0 software (Voorrips, 2002).

RESULTS ddRAD Sequencing and Identification of SNP Markers
The ddRADseq protocol was used to construct reduced representation libraries for the parents (Yuanza 9102 and Xuzhou 68-4) and their 188 RIL lines. Genomic DNA was double digested separately with restriction enzymes EcoR I (GAATTC) and Mse I (TTAA). Size of insert fragments in the libraries varied from 160 to 360 bp. The ddRAD-seq libraries were sequenced on the Illumina HiSeq4000 platform. In total, ∼361.03 Gb clean data (Q20>96%) were generated, containing 3719.77 million reads, with each read being ∼90 bp in length. Among these data, 161.39 million reads came from the parents, and ∼3558.38 million reads came from the libraries for the 188 RIL lines. The reads numbers of RIL individuals were 4.99-40.68 million range, with an average of 18.93 million reads (Figure 1).
After ustacks clustering analysis and SNP detection of all clean reads, 9,688 and 20,728 SNPs were identified from the reads of the specific locus EcoRI and MseI flanking genomic DNA fragments respectively. Among these, a total of 18,494 homozygous and polymorphic SNP markers were detected. The obtained stacks of each RIL lines were aligned against the catalog of all loci of two parents to determine their genotype at each locus. To avoid interference of paralog and repetitive sequence, the consensus sequence of the obtained SNPs were blasted against AA and BB genome of wild peanut and only the uniquely mapped reads were retained. Finally, 2,324 makers in the genome were identified for map construction.

Construction of High-Resolution Genetic Map and Comparative Analysis
The genotypes of the 188 RIL F 7 individuals were used to construct the genetic map using Joinmap 4.0 software. Overall, a genetic map was constructed, which contained 2,187 SNPs LGs, representing 20 chromosomes, with a map length of 1,566.1 cM (Figure 2 and Table 1). The SNP marker sequences on the genetic map were listed in Supplementary Table S1. The LGs ranged from 59.8 to 104.6 cM in length with an average of 78.30 cM, and 10 LGs contained over 100 markers. A02 is the shortest LG spanning 59.8 cM and contains only 40 markers, whereas A01 is the longest LG spanning 104.6 cM and contains 141 markers. The marker densities ranged from 0.37 to 2.22 cM/locus, resulting in an average distance of 0.72 cM between markers for the entire map (Figure 2 and Table 1). The Chi square analysis identified 987 SNPs (45.1%) with segregation distortion. All 20 LGs had good syntenic relationship with correspondence to AA genome and BB genome of diploid wild peanut (Figure 3).

Phenotypic Variation of Bacterial Wilt Resistance
The survival rate of parents and their RIL lines were investigated in disease nursery under four different environment conditions (2014HA, 2015HA, 2015NC, and 2016NC). The survival rate of the resistant parent Yuanza 9102, was significantly higher than that of the susceptible parent Xuzhou 68-4 (Figure 4 and Table 2). Frequency distribution showed that BW resistances had good consistency in four environment conditions ( Figure 5 and Table 2). Obvious phenotypic variations of resistance to BW were observed among RILs population under all four environment conditions. The survival rates of individual lines ranged from 0 to 100% and displayed two-peak distribution in four environment conditions (Figure 5). These results suggest that the resistance to BW is controlled by a major effect gene and a series of minor effect genes. ANOVA of BW trait across all four environment conditions indicated that genotype (G), environment (E) and genotype-environment interactions (G × E) had significant effects on the trait (Table 3). Moreover, the trait of BW resistance showed high broad-sense heritability of 81.72% (Table 3). These results suggest that genetic factor plays a major role in the expression of BW resistance.

Detection of QTLs for Bacterial Wilt Resistance
To identify QTLs for BW resistance, phenotypic data of four environment conditions were analyzed together with the genetic mapping data of 2,187 markers. A total of four significant QTLs were identified through single-environment QTL analysis. The four QTLs were all distributed on chromosome B02 and could explained 7.72-23.33% phenotypic variation of resistance to BW (Figures 6, 7 and Table 4). Among them, the stable QTL qBWB02.1 was repeatedly detected in three environment conditions of 2015NC, 2016NC, and 2015HA (Figures 6, 7 and Table 4). The major stable QTL (qBWB02.1) could explain phenotype variation of 12.17-23.33%, and the highest PVE Left three lines were Yuanza 9102, which is the BW-resistant parent. Right three lines were Xuzhou68-4, which is the BW-susceptible parent.
(23.33%) with LOD of 13.29 was in 2016NC environment. The qBWB02.1 was located at 0-3.5 cM map position on chromosome B02 with flanking markers AHMS173054 and AHES415390 (Figures 6, 7 and Table 4). This region could be traced to the genomic region of B02 of A. ipaensis. Corresponding position of 3.5 cM length on the genetic map was about 2.3 Mb in

Prediction of Candidate Resistant Genes on Chromosome B02
To identify the candidate genes for BW resistance in the stable major QTL loci (qBWB02.1), the target region (AHMS173054 to AHES415390) was examined for predicted genes according to A. ipaensis reference genome annotation database at Peanutbase 1 and the analysis of nucleotide-binding siteleucine-rich repeat (NBS-LRR)-encoding genes of A. ipaensis reference genome by Bertioli et al. (2016). The region of qBWB02.1 spanning 3.5 cM linkage interval corresponds to ∼2.3 Mb genomic region of B02 (physical position B02: 2,501,128-4,797,039 bp) of A. ipaensis. A total of 154 annotated genes were identified in this region. Of which, 21 genes (Table 5) were found to have disease resistance protein function. All these genes belonged to NBS-LRR genes, which was the largest class of resistance (R) genes. The 21 resistance proteins contained NBS domain and without Toll-interleukin-1 receptor (TIR) domain. Leucine-rich repeat (LRR) domain was also included in all proteins, except Araip.8GZ6F (Table 5). Our data suggests that the major QTL qBWB02.1, which was consistent and stable in three environment conditions, was a resistant genes-dense region. Therefore, these 21 genes harboring this genome segments provide potential disease resistant genes for BW resistance and warrant further investigation.

DISCUSSION
Cultivated peanut (A. hypogaea L.) is an allotetraploid (AABB, 2n = 4 × = 40) species with a large genome and low genetic diversity within the species. Because of these features, the construction of genetic linkage map in cultivated peanut has always been a formidable task. In recent years, many efforts have been made to construct genetic maps (Halward et al., 1993;Burow et al., 2001;Herselman et al., 2004;Foncéka et al., 2009;Varshney et al., 2009;Hong et al., 2010;Gautami et al., 2012;Qin et al., 2012;Shirasawa et al., 2012Shirasawa et al., , 2013Sujay et al., 2012;Wang H. et al., 2012;Seck et al., 2013). (2012) generated a map covering 2,166.4 cM with 1,114 loci using SSR and transponson markers. However, these maps were constructed using low throughput molecular markers with low density, thus unable to provide precise information on the QTLs controlling the traits of interest. In addition, the development of these polymorphic markers requires more screening work, cost for primer synthesis as well as resources and time. More efficient type of marker and a high amount of markers are crucial for improving the density and quality of linkage maps in peanut. SNPs are particularly attractive for genetic map construction, as they represent the most common type of DNA polymorphism in the genome and are amenable to high-through genotyping. Currently, there are a few studies which develop the SNP markers in peanut using the method of reduced genome sequencing, chips and resequencing (Nagy et al., 2012;Zhou et al., 2014;Shasidhar et al., 2017;Agarwal et al., 2018). Using ddRAD-seq method, Zhou et al. (2014) constructed the first high density SNP-based linkage map in tetraploid peanut consisting of 1,685 markers, 1267 loci, spanning a distance of 1,446.7 cM. Using resequencing method, Agarwal et al. (2018) developed the densest genetic map currently available for cultivated peanut with 8,869 SNPs and 2,156 mapped loci. In our SNP-based genetic map, 2,187 SNPs (1,962 loci) were identified using ddRAD-seq and were assigned to 20 LGs corresponding to 20 chromosomes of cultivated peanut. The total length of the genetic linkage map was 1,566.1 cM, with an average inter-loci distance of 0.72 cM. The genetic map constructed in this study is highly costeffective. On the other hand, the linkage map constructed in this study used a RIL F 7 generation as material. In general, the markers were more likely skewed segregation in high generation population of RIL population. This phenomenon probably related to gametophyte and/or zygotic natural selection and chromosomal rearrangements in many generations and artificial sampling involved in the construction of the RIL population. In fact, 987 (45.1%) markers of this genetic map showed skewed segregation. While the high generation population could improve the accuracy of bioinformatics analysis for SNP discovery. On the whole, the big loci number of the genetic map and its good synteny with AA and BB genomes provide a good foundation for important traits QTL mapping in peanut.
Based on this high-density genetic map, we detected four QTLs associated with BW resistance, including a major QTL (qBWB02.1) and three minor QTLs (qBWB02.2, qBWB02.3, qBWB02.4). qBWB02.1 was detected in both locations (Nanchong and Hongan) and 2 years (2015,2016) and was considered to be consistent and stable. qBWB02.3 and qBWB02.4 were detected in 2 years (2014,2015) and were considered to be consistent. Comparing the QTLs for BW resistance identified with previous studies (Jiang et al., 2007;Ren et al., 2008;Zhao et al., 2016), the QTLs located in our study were shown to be novel QTLs for BW resistance and the QTL interval in our study was significantly smaller than previously studies.
All four QTLs had positive additive effects, indicating these alleles were derived from resistant parent Yuanza 9102. From the pedigree of Yuanza 9102, we knew that the BW resistance was introgressed from diploid wild species of A. chacoense. A. chacoense is a good resistance resource of BW and some high BW-resistance accessions were obtained from this species (Chen et al., 2008). On the other hand, some different resistant alleles in QTLs detected by Zhao were derived from Yueyou 92 (Zhao et al., 2016). The source of resistance of Yueyou 92 can be traced back to Xiekangqing, a famous resistant cultivated peanut in China. In addition, previous study showed that the resistant gene of Yuanza 9102 exhibited dominance or partial dominance effects , while the resistant gene of Yueyou 92 was recessive (Zhao et al., 2016). The different QTL results may come from the different materials used.
In the present study, 21 disease resistance genes (Table 5) in the region of major QTL qBWB02.1 were identified. We considered these genes as putative candidate genes for BW resistance. All proteins encoded by the 21 genes contained the necessary domains (NBS and LRR) of NBS-LRR class. The domain can induce a programmed cell death or hypersensitive response (HR) in the presence of a specific pathogen elicitor (Bonardi et al., 2011). Although several QTLs related to resistance to BW have been identified in tomato (Wang et al., 2000;Carmeille et al., 2006), tobacco (Qian et al., 2013), Arabidopsis thaliana (Godiard et al., 2003), Medicago truncatula (Ben et al., 2013), very few resistance genes have been identified except for ERECTA gene and RRS1-R gene of A. thaliana. RRS1-R gene in A. thaliana, which is involved in monogenic resistance to BW, contains typical TIR-NB-LRR domain. It recognizes avirulence gene PopP2 and traffics to the nucleus through nuclear localization signal (Deslandes et al., 2002(Deslandes et al., , 2003. Another R gene ERECTA in A. thaliana, as a  quantitative resistance locus, involves in polygenic resistance to BW. In addition, overexpression of peanut NBS-LRR gene AhRRS5 in tobacco significantly enhanced the resistance to BW (Zhang et al., 2017). The QTL qBWB02.1 identified in this study was repeatedly detected and explained 12.17-23.33% phenotype variation of BW resistance in three environment conditions. Notably, the cluster of the resistant genes in this region implied that this genome region has potential gene resources for BW disease resistance and worthy further study.

CONCLUSION
We used ddRAD-seq technology for large-scale identification of SNPs that were subsequently used for high-throughput genotyping and construction of a high quality genetic map of cultivated peanut. Combining the phenotype data of four environment conditions, we identified four QTLs located on B02 for BW resistance. The major stable QTL region (qBWB02.1) covered 2.3 Mb physical distance containing 21 NBS-LRR genes. These genes were considered as putative candidate resistant genes for BW. This study provides a major stable QTL for fine-mapping and contributions to peanut breeding for BW resistance through MAS.

DATA AVAILABILITY
The datasets generated for this study can be found at https://www.ncbi.nlm.nih.gov/sra/SRP158656.