QTL mapping of web blotch resistance in peanut by high-throughput genome-wide sequencing

Web blotch is one of the most important foliar diseases worldwide in peanut (Arachis hypogaea L.). The identification of quantitative trait loci (QTLs) for peanut web blotch resistance represents the basis for gene mining and the application of molecular breeding technologies. In this study, a peanut recombinant inbred line (RIL) population was used to map QTLs for web blotch resistance based on high-throughput genome-wide sequencing. Frequency distributions of disease grade and disease index in five environments indicated wide phenotypic variations in response to web blotch among RILs. A high-density genetic map was constructed, containing 3634 bin markers distributed on 20 peanut linkage groups (LGs) with an average genetic distance of 0.5 cM. In total, eight QTLs were detected for peanut web blotch resistance in at least two environments, explaining from 2.8 to 15.1% of phenotypic variance. Two major QTLs qWBRA04 and qWBRA14 were detected in all five environments and were linked to 40 candidate genes encoding nucleotide-binding site leucine-rich repeat (NBS-LRR) or other proteins related to disease resistances. The results of this study provide a basis for breeding peanut cultivars with web blotch resistance.

regarded as one of the most urgent issues to be addressed in some peanut growing areas.
The fungus P. araehidieola pathogens normally overwinter within the crop residue and plants of all ages are susceptible [6]. Cloudy weather with frequent rains and temperatures (15-30°C) favor fungal activity [14]. Spore germination and cuticle penetration occurred on leaflets of peanut and the symptoms appeared seven to 9 days after inoculation. The initial symptoms of small, irregular, brown to reddish brown lesions along the midrib of both young and old leaves were observed and the severe defoliation was observed in the infected fields of susceptible cultivars [15]. Following penetration of the cuticle, fungal hyphae grew just under the cuticle for a few mm beyond the point of penetration and ceased growth, which indicated that the defense responses of P. araehidieola infection is hypersensitive-type reaction [6].
Although the web blotch disease has considerable economic importance in peanut farming systems, previous researches on this disease are relatively few compared with those on other main peanut foliar diseases. Zhang et al. (2019) presented the draft genome sequence of a Phoma araehidieola isolate named Wb2 and indicated that the draft genome of Wb2 was about 34.11 Mb and contained 37,330 open reading frames (ORFs), with G + C content 49.23% [16]. Smith et al. (1979) reported that Virginia and runner market-type peanut cultivars are more resistant to web blotch than the Spanish markettype cultivars [17]. Zhang et al. (2011) showed that resistance to web blotch was controlled by three major genes and several minor genes, and identified one quantitative trait locus (QTL) located on linkage group (LG) 7 [18]. At present, researches on peanut web blotch mainly focus on the classification status of the pathogen asexual generation, the pathogen molecular biology, chemical control and epidemic rules [6-9, 14, 15]. A similar disease called net blotch, caused by Pyrenophora teres is also one of the major diseases in barley [19]. Many studies have identified a large region of chromosome 6H responsible for resistance to net blotch in barley [20][21][22][23] and a region in chromosome 3H harbors two predicted genes from the NBS-LRR gene family [24].
With advances in genomic sequencing technologies and the availability of diploid and tetraploid genome assemblies in Arachis species [1,4], high-throughput genome-wide sequencing has become a primary strategy in peanut to identify single nucleotide polymorphisms (SNP) markers linked to resistance genes and QTLs. For example, four bacterial wilt resistance QTLs were identified on chromosome B02 using a RIL population and 2187 SNP markers [25]. A major QTL for resistance to late leaf spot, on chromosome B05, and two major QTLs for resistance to early leaf spot, on chromosomes A03 and B04, were mapped using a high-density genetic map comprising 2753 SNP markers and a F 9 RIL population of 192 individual lines [26]. Additionally, 62 QTLs for 14 yield-related traits were detected on 12 chromosomes across three environments using a high-density genetic map including 2636 recombination bin markers and a F 6 RIL population of 242 lines [27].
In the present study, high-throughput genome-wide sequencing technology was used to obtain SNP markers, and a SNP-based genetic linkage map was constructed to identify the QTLs for resistance to peanut web blotch. The results of this study will help to better understand the mechanisms of interaction between P. araehidieola and peanut and to develop resistant varieties.

Evaluation of web blotch resistance
Wide phenotypic variations in response to peanut web blotch were observed among RILs under all the five test environments in this study and the two parental varieties of the RIL population also showed significant difference in web blotch resistance (Fig. 1). The distributions of the disease scale recorded during the years 2007 and 2008 were shown in Fig. 1a, whereas the distributions of the disease index recorded during the years 2012 and 2018 (both in the field and indoor) were shown in Fig. 1b The frequency distribution of disease index in the year 2012 was almost a normal distribution (Fig. 1b). The frequency distributions under other conditions displayed skewed distribution (Fig. 1a, c and d), which had no impact on the following mapping results because only the random error term of the phenotypic data was required to follow the normal distribution [28].

Sequencing, SNP and bin markers discovery
A whole-genome resequencing strategy was applied to construct paired-end libraries for the parental lines and their 212 RIL progenies. The length of DNA fragments in the libraries was about 350 bp. Approximately, 490 Gb of clean data (Q20 > 96%) were produced, resulting from 6285 million reads, each with a length of 150 bp. In total, 600 million reads were generated for each of the two parents, whereas reads generated for each of the RILs varied from 22.37 to 24.83 million (Supplementary Table S1). Coverage rate, mapped reads rate, sequence depth and other results from alignment to the reference genome were shown in Supplementary Table S1. In particular, the coverage rate associated with the two parents Zheng8903 and Yuhua4 were 98.51 and 99.05%, respectively, whereas it ranged from 53 to 63.63% in the RIL population (Supplementary Table S1). The sequencing depth was 35.23× for both parents, whereas it ranged from 1.31× to 1.46× for the RIL population. Originally, 636,831 SNPs were called from the 214 samples using the GATK protocol. Then, 556,615 SNPs were retained after filtering for low quality loci in the two parents, due to missing values, heterozygosis, depth < 10 and GQ < 20. Finally, 138,039 SNPs which were homozygous and polymorphic between two parents were used for further analyses.

Construction of physical recombination maps and high density genetic linkage map
To avoid errors caused by low coverage associated with RIL sequencing, a sliding window with 15 consecutive SNPs was used to find the more accurate recombined breakpoints. The physical recombination map of 212 RILs was constructed based on the recombination map of each progeny (Supplementary Figure S1). After that, all chromosomes of the 212 RILs were aligned and compared for the minimal of 100-kb intervals. As a result, a total of 3634 bin markers for the 212 lines were obtained in this way, and the genotypes and physical locations of the bins were given in Supplementary Table S2.
The obtained 3634 bin markers were used to construct a genetic linkage map by the software JoinMap®v5.0 [29]. Twenty linkage groups were generated and assigned to the 20 chromosomes of the cultivated peanut according to the physical positions. The total genome length was 1817.91 cM and the marker density across the 20 linkage groups ranged from 0.39 to 0.66 cM with an average of 0.50 cM ( Table 1). The LG16 had the lowest marker number (129) and the shortest genetic length (54.58 cM), while LG3 had the highest marker number (277) and the longest genetic length (135.61 cM) ( Table 1 and Fig. 2). More than 97.5% of the inter-markers distances were lower than 3 cM. The highest inter-marker distance (16.06 cM) was associated with LG6 (Table 1).
QTL mapping and candidate genes prediction for peanut web blotch resistance QTL mapping of peanut web blotch resistance was performed with MapQTL® v6.0 [30], using phenotypic data collected across five environments. Eight QTLs associated with peanut web blotch resistance, located in eight different LGs, were confirmed in at least two environments, explaining from 2.8 to 15.1% of phenotypic variation and displaying LOD values ranging from 1.32 to 7.45 (Table 2). Two QTLs (qWBRA04 and qWBRA14) located on LG04 and LG14 were significantly associated with resistance in all the five testing environments in this study ( Table 2, Fig. 3 and Fig. 4) and explained more than 10% of phenotypic variation, indicating they are probably the major QTLs with stable expression. Except for qWBRA13 and qWBRA05, which were detected in four and two environments respectively, the other four QTLs qWBRA03, qWBRA16, qWBRA17, qWBRA19 were detected in three environments ( Table 2). Absolute values displayed by the additive effect ranged from 0.11 to 8.29, and were negative for the QTLs on LG4, LG5, LG13, LG14, LG19 (indicating that the favorable alleles originated from the resistant parent Zheng8903), and positive for the QTLs on LG3, LG16 and LG17 (indicating that the favorable alleles originated from the susceptible parent Yuhua4).
To identify candidate genes for peanut web blotch resistance, coding sequences in the genomic region associated with the QTLs qWBRA04 and qWBRA14 were examined for predicted function, according to the Arachis hypogaea cv. Tifrunner reference genome annotation database [4]. A total of 40 candidate genes were identified with a putative role in disease resistance (Table 3). In detail, the region of qWBRA04, spanning a   such as a bZIP transcription factor and a WRKY transcription factor-like protein.
To further validate the candidate genes, the mutation type of the SNPs located in the candidate genes were analyzed (Table 4). There were 24 polymorphic SNPs between two parents located in seven genes, eight of which were in noncoding DNA regions and two of which were synonymous mutation, and the remained 14 SNPs on chromosome A03, A04, and A14 were non-synonymous. The 14 non-synonymous SNPs were located in four coding sequence regions, four of which were in the region of Arahy.LFE0TK, one in Arahy.Q7VTCQ, three in Table 2 Quantitative trait loci (QTLs) for peanut web blotch resistance identified in this study QTL LG  Arahy.1RZ0PJ and six in Arahy.MR7539. All of the four genes were NBS-LRR, and three of them contained TIR domain except that Arahy.Q7VTCQ contained CC domain.

KASP design and validation
To realize molecular marker assisted breeding for web blotch resistant peanut varieties, the 24 polymorphic SNPs between two parents (  Table S3).

Discussion
A high density genetic linkage map with 3634 bin markers was constructed based on high-throughput whole-genome sequencing and a sliding window strategy. The order of the markers on the genetic map was overall consistent with the physical order according to the peanut genome assembly, except for some translocation occurring between LG3 and LG13 and between LG6 and LG16. A total of 277 bin markers were distributed on LG3, while the physical positions of the last 70 markers were on chromosome A13, and the physical positions of the last 47 markers of LG13 were on chromosome A03. Moreover, 43 bin markers from chromosome A16 were inserted in the LG6. The reason for such seeming translocation might be that there were some assembly errors in the reference genome which were illustrated in Peanutbase (https://www.peanutbase. org/data/public/Arachis_hypogaea/Tifrunner.gnm1.KYV3/) [31]. The genetic orders of the markers on LG1, LG4, LG5, LG10, LG12 and LG18 were fully consistent with their physical orders and a few markers on the other LGs were inversed with the adjacent markers. In total, eight novel QTLs distributed on eight chromosomes were identified for peanut web blotch resistance, which were confirmed in at least two of the experimental trials. Three of them were major QTLs, as they were associated with phenotypic variance explained (PVE) > 10%, whereas the other five were minor QTLs, in accordance with the predicted results of Zhang et al. (2011). Except for qWBRA03, which was detected only in naturally conditions, the other seven QTLs were detected both in the natural and artificial inoculation conditions. The two major QTLs qWBRA04 and qWBRA14 were stably detected across all the five testing environments in this study, and thus might be of great potential in breeding resistant peanut varieties.
In the present study, we identified 55 candidate disease resistant genes in the target region of eight QTLs  (Tables 3 and 5), of which 40 were linked with the two major and stable QTLs qWBRA04 and qWBRA14 (Table  3), while the other 15 were associated with the other six QTL intervals ( Table 5). The 40 candidate genes covered by the intervals of qWBRA04 and qWBRA14 included 21 TIR-NBS-LRR and 3 CC-NBS-LRR, which are the two well-known R gene types [32]. Also, there were two candidate genes encoded LRR and NB-ARC (nucleotidebinding domain shared with APAF-1, various R-proteins and CED-4) domain, which also encoded resistance genes [33]. Among the remained 15 genes, seven of them function in the downstream pathways of resistant signaling, the other eight encoded proteins contain disease resistance response related domains but could not assigned to the well-known R-gene types (Table 3). A total of 26 NBS-LRR genes were identified to be related to peanut web blotch resistance in this study. NBS-LRR is the biggest category of R genes [34] and has been identified at the genome-wide level in Arachis [35]. It had been found that NBS-LRRs were involved in response to late leaf spot, tomato spotted wilt virus, and bacterial wilt in A. duranensis, A. ipaensis, and A. hypogaea [36]. It could be concluded that NBS-LRR was also involved in the resistance to peanut web blotch, but the regulatory mechanism in the process of disease resistance needs to be further studied.
The results of three validated KASP markers indicated that the gene Arahy.Q7VTCQ (CC-NBS-LRR) might be one of the resistant genes for peanut web blotch in a great possibility. The reason of the inconsistence between the phenotypes and genotypes of three test lines may be that the three KASP markers were not completely linked with peanut web blotch resistance. Therefore, further study will be needed to design closer linked markers to be employed in molecular marker assisted breeding (MAS).

Conclusion
In this study, eight QTLs for peanut web blotch resistance were detected and two major QTLs qWBRA04 and qWBRA04 were linked to 40 candidate genes encoding NBS-LRR or other proteins related to disease resistance, which may shed some insights on understanding web blotch resistance and facilitate the development of resistant peanut cultivars.

Plant materials
A RIL population consisting of 212 F 12 lines derived from the parental cross combination between lines Zheng8903 and cultivar Yuhua4 was used in this study. The female parent Zheng8903 with the pedigree '79-266//71-31/Chico' is a breeding line showing high level of resistance to peanut web blotch [37]. The male parent Yuhua4 is a variety released in 1991 by the Henan Academy of Agricultural Sciences and we obtained the seeds from our own inventory [38], showing high level of susceptibility to web blotch [37]. All the plant materials including RILs and its parents mentioned above were developed and preserved in the corresponding author's lab.

Experimental trials and phenotyping
All the experimental trials to evaluate response to peanut web blotch were conducted at the research station of the  [39]. For field inoculation in the year 2018, plots were sprayed with an inoculum of 1.6 × 10 − 3 g/ml at the flowering stage [40] and disease evaluation was carried out 20 days after inoculation according to the 0-9 scale described in Yu (2011) [41]. For indoor inoculation, five plants were inoculated for each line with two replications at the 6 leaves stage. The inoculum concentration was 2 × 10 6 conidia/ml and the preparation of conidia suspension described by Zhang (2019) [37] was followed. Two weeks after inoculation, 12 inoculated leaves at the main stem were collected and the lesion area for each was scanned by the Leaf Area Meter (Wanshen LA-S).

Sequencing and genotyping
Genomic DNA of the parents and 212 RILs was extracted from young leaf tissues using the Plant genome DNA extraction kit (TIANGEN) and randomly sheared by sonication. The DNA fragments with the length of 300 bp were recovered by electrophoresis. After ligating DNA fragments with adapters, libraries were paired-end sequenced using the Illumina Hi-seq platform with read length of 150 bp. After trimming adapters and low quality reads, clean data was used for aligning to the reference genome, allowing SNP identification and genotyping. In detail, the assembly of Arachis hypogaea cv. Tifrunner was used as the reference genome [4]. The aln command in the software bwa-0.7.10 was used to align clean data to the reference genome, and unique reads were used for subsequent SNP variation detection by the software Linkage map construction and QTL mapping As the accuracy of call at single SNP loci was low, due to the low coverage chosen for RIL sequencing, a sliding window approach was applied to evaluate a group of consecutive SNPs for genotyping [42]. The genotypes of all chromosomes of the 212 RILs were aligned and compared for the minimal of 100-kb intervals and the adjacent 100-kb intervals with the same genotype across the entire RIL population were recognized as a single recombination bin [42]. From bin markers, linkage groups were constructed using JoinMap v5.0 [29], selecting LOD scores from 2 to 10 to identify groups and the regression algorithm to perform ordering within each LG. The final linkage map was drawn using the R package LinkageMapView. QTL mapping was performed using MapQTL v6.0 [30], by selecting multiple QTL mapping (MQM) to detect potential QTLs with the LOD threshold of 2.5 in at least one environment.