QTL mapping and identification of genes associated with the resistance to Acanthoscelides obtectus in cultivated common bean using a high-density genetic linkage map

Common bean (Phaseolus vulgaris L.) is an important agricultural product with large nutritional value, and the insect pest Acanthoscelides obtectus (Say) seriously affects its product quality and commodity quality during storage. Few researches on genes of bruchid resistance have investigated in common bean cultivars. In this study, a bruchid-resistant cultivar black kidney bean and a highly susceptible accession Longyundou3 from different gene banks were crossed to construct a recombinant inbred line population. The genetic analysis indicated a quantitative inheritance of the bruchid resistance trait controlled by polygenes. A high-density genetic map of a total map distance of 1283.68 cM with an average interval of 0.61 cM between each marker was constructed using an F6 population of 157 recombinant inbred lines. The map has 3106 bin markers, containing 2,234,769 SNPs. Using the high-density genetic map, a new quantitative trait locus for the resistance to Acanthoscelides obtectus was identified on chromosome 6. New molecular markers based on the candidate region were developed, and this locus was further delimited to an interval of 122.3 kb between SSR markers I6–4 and I6–16 using an F2 population. This region comprised five genes. Phvul.006G003700, which encodes a bifunctional inhibitor, may be a potential candidate gene for bruchid resistance. Sequencing analysis of candidate gene identified a 5 bp insertion-deletion in promoter of gene Phvul.006G003700 between two parents. Expression analysis of candidate gene revealed that the expression level of Phvul.006G003700 in bruchid-resistant parent was markedly higher than that in bruchid-susceptible parent both in dry seeds and leaves. A high-density genetic linkage map was constructed utilizing whole-genome resequencing and one new QTL for bruchid resistance was identified on chromosome 6 in common bean cultivar. Phvul.006G003700 (encoding a bifunctional inhibitor) may be a potential candidate gene. These results may form the basis for further research to reveal the bruchid resistance molecular mechanism of common bean.


Background
Common bean (Phaseolus vulgaris L.) is one of the most widely adapted legumes in the developing world and the most important food legume for direct human consumption in the world [1,2]. This autodiploid (2n = 2x = 22) species has 11 chromosomes and a genome size of about 537 Mb (Phaseolus vulgaris v2.1, DOE-JGI and USDA-NIFA, http:// phyto zome. jgi. doe. gov/). Common bean is an annual and self-pollinating plant, and its seeds are rich in protein, amino acids, flavonoids, alkaloids and terpenoids, which are necessary for human life [3][4][5][6]. Common bean originated from Mesoamerica, and China is considered as a secondary centre of common bean diversity and has a long history of cultivation [7,8]. Common bean is an important source of plant protein in the human dietary structure and is one of the main export products of China [9].
During their storage, legume seeds are often damaged by legume weevils, which lead to a serious decrease in yield. Bean weevils [Acanthoscelides obtectus (Say)] and Mexican bean weevils [Zabrotes subfasciatus (Boheman)] are the main insect pests of common bean during storage [10,11]. Acanthoscelides obtectus (Say) originated from South America and Central America and has been introduced into dozens of countries and regions, such as Afghanistan, Japan, Myanmar and China, through international trade and introduction channels [12]. The larvae of A. obtectus can enter the bean seeds for feeding. The adults fly to the fields or warehouses and continue to damage common bean, causing great losses [13,14]. The damage and weight loss rates of common bean in the field are 30% ~ 38% and 6% ~ 8%, respectively, and storage increases these values to 74% and 64%, respectively, which indicates that the damage created by bruchids is more serious during the storage period [15][16][17]. At present, the commonly used physical methods, including freezing [18][19][20], boiling and seed dressing treatment [21] as well as chemical methods of fumigation with chemical agents such as phosphine [22] and ethyl formate [23], and biological control [24,25] are applied to control seeds damage by bruchids, but the negative impacts on the seeds themselves and the environment are substantial.
At present, researches on bruchid resistance genes cloning are lagging behind both at home and abroad, and in particular, few studies have investigated on resistance genes in common bean cultivars using molecular markers. Blair et al. [44] used SSR molecular markers to identify a locus linked to arcelin 1 allele from the wild common bean accession RAZ106 and found that this locus was located at the APA locus on Pv04. Kamfwa et al. [45] mapped two QTLs on chromosome 4 (one associated with the APA locus) and one QTL on chromosome 6 using an A. obtectus resistant breeding line introgressed from the wild tepary bean (contained Arcelin gene). Lectin can be encoded by a single locus or by several loci, for example, APA locus located on Pv04, Lec-2 and Lec-3, and the lec genes mapping on Pv07 [46,47].
In this study, we found that the black kidney bean, a cultivar of common bean, has high resistance to Acanthoscelides obtectus (Say). We then constructed a highdensity genetic linkage map using a recombinant inbred line population obtained from a cross between the bruchid-susceptible cultivar Longyundou3 and black kidney bean and performed gene mapping of bruchid resistance. This study mapped and analyzed the resistance genes for bean weevils in cultivated common bean, which may lay a foundation for further study on the mechanism of bruchid resistance in common bean and provide a possibility for improved seed preservation.

Bruchid resistance segregation and genetic analysis
The black kidney bean (BKB) exhibited resistance to Acanthoscelides obtectus (Say), with a mean percentage of damaged seeds of 28.3%. In contrast, Longyundou3 (LYD3) showed high susceptibility, and all the seeds were damaged by bruchids (percentage of damaged seeds equal to 100%) (Fig. 1a, b).
The resistance of F 2:3 population and RILs was identified. Two traits, namely, the percentage of damaged seeds (PDS) and the number of perforations (NP), were investigated. The PDS among the F 2:3 families ranged from 0 to 100%, and the NP ranged from 0 to 220, showing significant resistance segregation (Fig. 1c). The PDS in the RILs ranged from 0 to 100%, and the NP ranged from 0 to 280. However, the minimum damage rate of seeds was lower than that of resistant parent, which indicated that resistant trait had a positive superparent advantage. The distribution of bruchid resistance was continuous in the F 2:3 population, and the ratio of susceptibility to resistance was not consistent with Mendelian inheritance. The kurtosis and skewness values were − 0.84 and 0.11, respectively, with the absolute values lower than one, which was approximately normal, and these results indicated the quantitative inheritance of the resistance trait controlled by polygenes.

High-density genetic linkage map construction
Based on the technical characteristics of whole genome resequencing, the 157 RILs and two parents of common bean were resequenced. Raw reads were obtained by CASAVA base calling and conversion from the raw image data acquired by high-throughput sequencing. After filtering, 80,029,096 (17.12× coverage) and 71,913,622 (16.95× coverage) clean reads were generated from BKB (resistant parent) and LYD3 (susceptible parent), Fig. 1 Phenotypic characteristics and frequency distribution of PDS and NP in the F 2:3 families. The F 2:3 families were derived from a cross between LYD3 (bruchid-susceptible parent) and BKB (bruchid-resistant parent). a Phenotypic characteristics of the two parents before infestation with bean weevils. b Phenotypic characteristics of the two parents after infestation with bean weevils. The large-yellow seeds on the left are the susceptible parent LYD3, and the small-black seeds on the right are from the resistant parent BKB. c Frequency distribution of the PDS and NP in the F 2:3 families from a cross of LYD3 and BKB respectively. The clean reads from 157 RILs ranged from 39,686,552 to 105,924,272, with an average sequencing depth of 14.38× coverage. By comparing the sequencing data of BKB, LYD3 and 157 RILs with the reference genome, a total of 300,003 SNPs among the parents were screened out using GATK 3.8 software [48].
The sliding window approach, with 15 SNPs in each sliding window, was used in this study. According to the proportion of SNP markers in the sliding window from the male or female parents, the genotype of the segment in where the sliding window located was determined, and the location of the recombination breakpoint of the chromosome segment was determined. When the sliding window hit a homozygous/homozygous breakpoint, the genotype transiently changed from homozygous to heterozygous and then changed to the other homozygous genotype. A genotype remained unchanged until it hit a recombination breakpoint [49]. After all genotypes were called and the recombination breakpoints were determined, we identified a total of 48,689 recombinant breakpoints, with an average of 310 breakpoints per RIL (Fig. S1a).
By consolidating continuous non-recombination markers on the genome into a bin, a total of 4214 recombination bins were obtained for the 157 RILs. The average physical length of the recombination bins was 135.2 kb, including 712 SNPs. Each bin was used as a marker for genetic marker screening and genetic linkage map construction. Because segregation distortion could affect the results of map construction and QTL mapping, a ratio of 1:1 was used to perform segregation filtration with 4214 candidate markers, and 1108 biased markers were filtered out to yield 3106 effective markers. According to the genome information [50], the linkage groups were constructed using JoinMap 4.0 software [51] by controlling the LOD value between 2.0 and 10.0. A total of 11 linkage groups were generated with 3106 bin markers, containing 2,234,769 SNPs. Taking the linkage group as the unit, the genetic distance between the markers was calculated using the Kosambi algorithm. A high-density genetic map of a total map distance of 1283.68 cM, with an average interval of 0.61 cM between the bin markers, was constructed ( Table 1, Fig. 2). The number of effective bin markers in each chromosome ranged from 134 to 417. The length of each chromosome ranged from 79.23 cM to 166.73 cM, and the average distance between bin markers ranged from 0.43 cM to 0.80 cM. With the exception of biggest gap > 10-cM on LG05 and LG08, all other biggest gaps were < 8-cM. The conversion ratio between genetic distance and physical distance ranged from 1.67 to 2.79. LG08 was the longest with a map length of 166.73 cM, containing 417 bin markers with 297,249 SNPs, and had an average genetic distance of 0.71 cM.
LG03 was the second longest with a map length of 159.84 cM containing 289,659 SNPs. LG09 and LG01 Collinearity analysis was performed based on the positions of markers on the genome and the genetic map, and the results showed that the order of most marker positions on each linkage group was consistent with the chromosome (Fig. S1b), indicating good

QTL analysis and mapping for bruchid resistance
The genetic map was used to identify QTLs controlling bruchid resistance in common bean. After the identification of bruchid resistance, the PDS of the RILs ranged from 0 to 100%, and the NP ranged from 0 to 280. There was a significant positive correlation between the two phenotypes, with a correlation coefficient of 0.88 (Fig. S2). The ICIM-ADD method in QTL IciMapping 4.2 software was used to analyse PDS and NP. The analysis of the 157 RILs detected one QTL associated with PDS located on chromosome 6 between markers Bin1565 and Bin1566. The LOD value was 5.49, explaining 16.10% of the phenotypic variation, and the additive effect was − 0.11. One QTL related to NP, which was also located on chromosome 6 between markers Bin1565 and Bin1566, had the LOD value of 4.31 and additive effect of − 27.48, and explained 16.37% of the phenotypic variation ( Table 2; Fig. 3a; Fig. S3). This QTL for both PDS and NP between markers Bin1565 and Bin1566 was 0.23 cM, and located at the physical position between 1,250,000 and 1,500,000 bp on chromosome 6 ( Fig. 3b).  We named the QTL qAO6.1, and the beneficial alleles at this locus came from BKB, the resistant parent. The qAO6.1 might be a major locus for resistance to Acanthoscelides obtectus.
For narrowing down the region identified by QTL-seq, we developed 44 SSR markers and 5 InDel markers (Table  S1) on chromosome 6 in candidate region to detect polymorphisms between BKB and LYD3. A total 11 SSR markers and 3 InDel markers on chromosome 6 were polymorphic between the two parents. The regional linkage map for chromosome 6 were constructed with polymorphic DNA markers using F 2 mapping population consisted of 188 individual plants. We identified ten recombinants from the fine-mapping population using seven molecular markers. The qAO6.1 locus was delimited to an interval of approximately 122.3 kb between SSR markers I6-4 and I6-16 (Fig. 3c).

Annotation and candidate gene prediction for candidate region
According to the common bean reference genome sequence of G19833v1.0, the region between markers I6-4 and I6- 16

Sequence and expression analysis of candidate gene
The coding DNA sequences (CDS) of Phvul.006G003700 in BKB and LYD3 were determined. The amplification of the candidate gene indicated that Phvul.006G003700 had a length of 999 bp in cDNA in the resistant genotype BKB and encoded 332 amino acid residues. No difference was found between BKB and LYD3 (Fig. S5). The results from the search of the SMART database showed that the Phvul.006G003700 protein contained an AAI domain in the region from amino acids 55 to 329 (Fig. S6).
To further examine the transcriptional regulation of this gene, we analysed the upstream promoter sequence of the Phvul.006G003700 gene using the PlantCARE database. The BLAST searches revealed that the amplified sequence included essential cis-regulatory elements of the promoter such as the core promoter element TATA and the common element CAAT that were highly conserved across the investigated species. We obtained the sequence of promoter region for Phvul.006G003700 gene. A 5 bp difference on 800 bp upstream of the start codon in the promoter region of the Phvul.006G003700 gene was found between LYD3 and BKB, which resulted in the absence of essential cis-regulatory element TATA box in the susceptible material (Fig. 4). TATA box was known to be an essential transcription regulator. The absence of the TATA box may reduce the efficiency of transcription and thus affect the gene expression of Phvul.006G003700 in the susceptible parent LYD3.
The qRT-PCR results showed that the expression level of Phvul.006G003700 in BKB was markedly higher than LYD3 both in seeds and leaves (Fig. 5), which revealed that the gene expression level was higher in the bruchidresistant accession than in the bruchid-susceptible accession. In seeds, gene Phvul.006G003700 in BKB was about 37 times more expressed than it in LYD3, and about 3.5 times more expressed in leaves.
Through the results of association study between SNP markers and two traits of PDS and NP in the natural Fig. 4 Sequence differences in promoter region of gene Phvul.006G003700 between BKB (resistant) and LYD3 (susceptible). There was a 5 bp difference on 800 bp upstream of the start codon in the promoter region of gene Phvul.006G003700 between LYD3 and BKB, locating on the essential cis-regulatory element TATA box of the promoter, resulting in the absence of TATA box in the susceptible parent. The TATA box was marked in red box. The red asterisks indicate the different base sites between two sequences population, we found that three significant SNPs associated with PDS (Pv06_1376148, P = 5.85 × 10 − 10 ; Pv06_1376347, P = 9.97 × 10 − 10 ; Pv06_1376967, P = 1.11 × 10 − 9 ) were localized in the 3′-UTR region of the candidate gene Phvul.006G003700 (Fig. S7). These findings further support the likelihood that Phvul.006G003700 was the candidate gene for bruchid resistance.
To further clarify the relationship between the Phvul.006G003700 gene and homologous genes from other crops, a BLAST search against the NCBI and Uni-Prot databases was performed to find the homologous protein of Phvul.006G003700 in soybean (Glycine max), lupine (Lupinus albus), wheat (Triticum aestivum), maize (Zea mays), rice (Oryza glumipatula), barley (Hordeum vulgare subsp. vulgare), Medicago truncatula, Arabidopsis thaliana and other crops. The protein sequence encoded by Phvul.006G003700 was compared with the homologous protein sequences in other crops by phylogenetic tree analysis (Fig. S9). The results showed that the Phvul.006G003700 protein was closely related to the AAI protein in soybean, which suggested that it might have a similar function to AAI. The analysis of conserved motifs revealed three motifs with an e-value >1E-60 among Phvul.006G003700 and homologous proteins in other crops, and the first two conserved motifs were located within the AAI domain region (Fig. S10).

Discussion
During the storage of common bean, the damage caused by the common bean weevils leads to a decrease in stored grain [16,17]. Although some physical and chemical methods [23,24,52] can be used to prevent and control the damage, the treated seeds are difficult to eat or sell again, resulting in relatively large economic losses [20]. The screening of resistant materials and identification of resistance genes of common bean are needed to provide a green and environmental protection approach. Previous studies have found some common bean germplasms that are resistant to bean weevils. Most resistant germplasms were concentrated in wild bean, and few were concentrated in cultivars/landraces. Based on the preliminary screening and identification of resistance to A. obtectus among several hundred common bean germplasm resources, we identified one cultivated common bean variety with good resistance to weevils. It was named black kidney bean (BKB), and presented an average seeds infestation rate of less than 40%. Using the cross between BKB and LYD3, a cultivar completely susceptible to A. obtectus, we constructed a segregated population. Through the identification of resistance to A. obtectus among the F 2:3 population, we demonstrated that the trait of bruchid resistance was a quantitative trait.
The genetic linkage map is a powerful tool for gene mapping and cloning and is generally constructed using various molecular markers, such as SSR and SNP [53]. In our study, a high-throughput method for genotyping recombinant populations utilizing whole-genome resequencing data [54][55][56][57] was implemented, and a sliding window approach was designed to detect genome-wide SNPs for genotype calling and recombination breakpoint determination. This method is efficient and fast, and widely used in r crops, such as wheat and cucumber [49,58]. We constructed a high-density genetic linkage map including 3106 bin markers using 157 common bean RILs derived from a cross between BKB and LYD3.The map contained 2,234,769 SNPs, and the number of markers were more than those in previous studies [59][60][61]. Compared to using F 2 mapping population, using an RIL as the mapping population in our study can effectively reduce the influence of the dominant effect and reveal the additive effect of QTLs [62].
The researches on bruchid resistance genes in common bean are relatively few at present, and most of the genes associated with bruchid resistance exist in wild species, such as Arcelin1 to Arcelin8 [34,36,39,40,[63][64][65][66]. In this study, we detected one QTL associated with both two traits, PDS and NP. This QTL mapped to the region on chromosome 6 between markers Bin1565 and Bin1566 and had the additive effects of − 27.48 and − 0.11, explaining 16.37% and 16.10% of the phenotypic variation, respectively. The additive effects were negative, indicating the beneficial alleles at this locus came from BKB, the resistant parent. We delimited this QTL to an interval of approximately 122.3 kb between SSR markers I6-4 and I6-16, and it may be a major locus for resistance to Acanthoscelides obtectus. This QTL was same located on chromosome 6, but approximately a distance of 22.9 Mb upstream of the QTL reported by Kamfwa [45], indicating that the qAO6.1 identified in this study was a new locus related to bruchid resistance.
The Phvul.006G003700 encodes a bifunctional inhibitor/seed storage/LTP family protein belonged to the AAI-LTSS protein superfamily, which is known to play important roles in defending plants from insects and pathogens, in lipid transport between intracellular membranes and in nutrient storage. The results from the search of the SMART database (https:// smart. embl. de) [67] showed that the Phvul.006G003700 protein contained an AAI domain. Previous studies have been reported that a double-headed inhibitor (subtilis protease/α-amylase inhibitor and trypsin/α-amylase inhibitor) in Ragi (finger millet, E. coracana Gaertneri) could simultaneously bind a molecular protease (subtilis protease and trypsin) and a molecular α-amylase [68][69][70]. The α-amylase was a kind of important digestive enzymes in the saliva and gut of insects [71,72], and proteinase inhibitors can combine with digestive enzymes to prevent insects from eating [73]. For example, proteinaceous inhibitor from finger millet could inhibit α-amylase of insect pests, and the highest inhibition percentage was found against α-amylase of Callosobruchus chinensis [74,75]. The α-amylase inhibitor is an important biochemical substance in common bean, and plays an important role in insect resistance. The previous studies demonstrated that the seeds of adzuki bean with a transferred α-amylase inhibitor gene could accumulate a high amount of α-amylase inhibitor and had a high resistance to bruchid [27]. Similarly, one transgenic chickpea and four cowpea lines expressing αAI-1 (an α-amylase inhibitor from common bean) were highly resistant to both bruchid pest species Callosobruchus chinensis and Callosobruchus maculatus [76]. Franco et al. [77,78] confirmed plant α-amylase inhibitors can inhibit the digestive enzymes from the bean weevil Acanthoscelides obtectus. According to the analysis of candidate gene sequences, we found that there was a 5 bp difference in the promoter region of Phvul.006G003700 between the two parents, which was just located in the TATA box region. TATA box is an essential transcription regulator and may affect gene expression level. The qRT-PCR analysis of gene Phvul.006G003700 showed that the expression level of Phvul.006G003700 in bruchid-resistant parent was much higher than that in bruchid-susceptible parent in dry seeds, which corresponded to the presence of TATA box in resistant parent and the absence of TATA box in susceptible parent due to 5 bp insertion. Basi´ studies [79] showed that the stepwise deletions of the TATA box of nmt1 gene highly expressed in Schizosaccharomyces pombe could result in the decrease in transcription efficiency and level of expression downregulated. Therefore, we speculated that the high level of Phvul.006G003700 gene expression in BKB could lead to the high content of α-amylase inhibitor, thus enhancing the inhibition to α-amylase of Acanthoscelides obtectus, and producing anti-insect effect. Further studies for the functional validation of this candidate gene are required, for example using complementation or gene knockout analyses, and physiological and biochemical researches.

Conclusions
We constructed a high-density genetic map by wholegenome resequencing using a RIL population. The map has 3106 bin markers, containing 2,234,769 SNPs. A new major QTL for resistance to A. obtectus on chromosome 6 was identified from a common bean cultivar based on the constructed genetic map. Phvul.006G003700 encodes a bifunctional α-amylase/protease-inhibited protein and is likely to be a candidate gene responsible for bruchid resistance. This research was very useful for screening and identifying bruchid resistance resources in common bean, and further research on the excavation of bruchidresistant genes.

Plant materials
The black kidney bean (BKB), a Chinese cultivar in the Mesoamerican gene pool of common bean that has a determinate bush plant and small black seeds (100-seed weight of 17.05 g ± 1.38), displays bruchid resistance. Longyundou3 (LYD3), a Chinese cultivar originating from the Andean gene pool, is an indeterminate viny plant with large yellow seeds (100-seed weight of 59.23 g ± 4.00) and susceptible to Acanthoscelides obtectus (Say). A genetic mapping population was created from a cross between LYD3 and BKB, and the F 6 population, which included 157 recombinant inbred lines (RILs), was obtained after continuous self-crossing. All seeds production was carried out at the experimental field of Heilongjiang Academy of Agricultural Sciences at the Harbin and Changping experimental sites of the Chinese Academy of Agricultural Sciences in Beijing, China.

Evaluation of bruchid resistance
The improved indoor artificial pest inoculation method was used to identify the resistance to Acanthoscelides obtectus. Acanthoscelides obtectus was collected from local red kidney beans infested with partially feathered bean weevils located in several cities and counties with severe bruchid damage in Yunnan Province. Forty healthy dry seeds were randomly selected from harvested seeds of each F 2 individual and F 6 line. Trials were replicated two times, with 20 seeds per replicate. Fifty healthy dry seeds of resistant parent and susceptible parent were used as control. Each material was placed into a 9-cm-diameter petri dish and randomly placed on the shelves. We covered all the shelves with a fine net, and put thousands of the collected feathered adults of A. obtectus (male and female) into the net to mate freely, so as to ensure that all the seeds were infested under the same condition. The room temperature was maintained at 18 ~ 25 °C. The humidity was maintained at approximately 70%, and timely ventilation was performed. After 100 days, the damage degree of each material was investigated. The percentage of damaged seeds (PDS) and the number of perforations (NP) were calculated. The percentage of damaged seeds was the percentage of infested seeds in the total number of identified seeds, and the number of perforations was the total number of infested holes of all identified seeds in each material.

Genetic map construction of the RIL population
BKB, LYD3 and 157 RILs from a cross between BKB and LYD3 were used for resequencing. The DNA samples of sufficient quality were randomly interrupted into 350bp fragments using a Bioruptor Pico ultrasonic crushing machine, and a genomic DNA library construction kit was used for library preparation, including DNA fragment end flatting, adding an ' A' to 3′ end, connecting an adaptor with DNA end and PCR amplification. Quality control was performed with the constructed libraries. The libraries of sufficient quality were sequenced on an Illumina HiSeq2500 for high-throughput sequencing, generally using PE150 sequencing mode, that is, double-end sequencing with 150 bp measured at each end. The raw image data files obtained by high-throughput sequencing were analysed by CASAVA base calling (http:// suppo rt. illum ina. com/ seque ncing/ seque ncing_ softw are/ casava. ilmn) [80] and converted into original sequencing sequences, which were called raw data or raw reads. After filtering, clean and effective reads were obtained. The clean reads of the two parents and 157 RILs were compared with the reference genome G19833v1.0 [50] using BWA software [81]. The single nucleotide polymorphisms (SNPs) of each sample were detected using GATK 3.8 software [48], and tens of thousands of SNP markers were developed.
The SNPs were filtered based on the missing percentage ≤ 20% and MAF ≥ 5%. We combined the SNP markers that were continuously unrecombined in the genome into a bin. First, the sliding window method was used to construct a recombination map. Each sliding window contained 15 SNPs. According to the proportion of SNP markers in the sliding window from the male or female parents, the genotype of the segment in where the sliding window located was determined, and the location of the recombination breakpoint of the chromosome fragment was identified. A recombination map was then drawn. Second, a bin map was constructed based on the recombination map of all the families in the group. Each bin was used as a marker for genetic marker screening and genetic linkage map construction. Taking the linkage group as a unit, JoinMap 4.0 software [51] was used to construct the genetic linkage map, and the Kosambi algorithm was used to calculate the genetic distance between markers.

QTL mapping for bruchid resistance
Quantitative trait loci (QTLs) were identified by inclusive composite interval mapping (ICIM) using QTL Ici-Mapping4.2 software [82]. A 1-cM scan window was employed, and the likelihood ratio statistic was computed every 2 cM. The LOD (log of odds) values and PVE were determined based on likelihood ratio tests under a hypothesis allowing both additive and dominance effects. A 1000 permutation test was used, and QTLs with LOD values of 3.0 and higher were called for candidate loci.

DNA markers development
The candidate region for bruchid resistance was verified using the analysis of SSR (simple sequence repeats) markers and InDel (insertion-deletion) markers. The SSR markers and InDel markers were designed in primer5 (Premier Biosoft International, USA). According to the location of candidate intervals in the reference genome G19833v1.0 [50], forty-four new SSR markers and 5 InDel markers on chromosome 6 were developed to screen for polymorphisms between BKB and LYD3 (Table S1). Polymerase chain reaction (PCR) was performed in 10 μL reactions containing 20-25 ng of template DNA, 0.2 μmol L − 1 forward and reverse primers each 0.5 μL, 5 μL of 2 × Taq DNA polymerase Mix and suitable nuclease-free water in a T100 Thermal Cycler (Bio-Rad Research, USA). The amplification profile was followed by 95 °C for 5 min, 30 cycles of 95 °C for 30 s, 50 °C-60 °C for 30s, 72 °C for 1 min/1 kb, and final extension for 10 min at 72 °C. The amplified products were electrophoresed on 8% nondenaturing polyacrylamide gels, and the bands were visualized by using silver nitrate staining.

Association study of candidate region
One natural population (a total of 628 common bean accessions representing broad genetic diversity, including 504 accessions in China and 124 elsewhere) was used in the association study for candidate region (Table S2). Genotypic data were generated from the whole genome resequencing project of 683 accessions as previously described by Wu [1], and nucleotide variations were further filtered based on the missing rates (≤ 10%) and a minor allele frequency of ≥ 0.05. Two traits of 628 common bean accessions, PDS and NP, were assessed in this study. The general linear model (GLM) [83] was used to detect the association between markers and phenotypes using Tassel 5.0. The P-value threshold was 1 × 10 − 6 [2], and the markers that exceeded the threshold (1 × 10 − 6 ) were considered to be associated with PDS or NP.

Annotation of candidate genes
The genes of candidate regions were annotated using the annotation file of common bean accession G19833v1.0 [50]. The homologous genes of common bean candidate genes were searched in Arabidopsis thaliana and other crops according to the NCBI (https:// www. ncbi. nlm. nih. gov/) and UniProt (https:// www. unipr ot. org/) databases.

Sequence amplification of candidate genes for bruchid resistance
Based on the results from QTL mapping, gene-specific primer pairs (Table S4) were developed and used to amplify the genomic DNA (gDNA), complementary DNA (cDNA) and promoter sequences of the candidate genes in BKB and LYD3. PCR amplification was performed in 50 μL volumes containing 100 ng of template DNA, 1 μL of each forward and reverse primer (10 μmol L − 1 ), 25 μL of 2 × Trans Start ® KD Plus Super Mix (−dye) and the suitable nuclease-free water in a T100 Thermal Cycler (Bio-Rad Research, USA). The thermal cycling programme was as follows: 94 °C for 5 min, 35 cycles of 94 °C for 30 s, 53-60 °C for 30 s, and 68 °C for 1 kb/min, and final extension at 72 °C for 8 min. The amplified products were electrophoresed on a 1.0% agarose gel. The target bands were cut and purified using a TIANgel Midi Purification Kit, then connected with the pEASY ® -Blunt Cloning vector and transferred into the newly thawed Trans1-T1 Phage Resistant Chemically Competent Cell. The connected product was coated on an LB plate containing IPTG, X-gal and ampicillin antibiotics and cultured overnight. White colonies were selected for detection and sequencing.

Bioinformatics analysis
Sequence alignment was performed by MEGA 6.0 software [84]. The domains of the proteins encoded by the candidate genes were identified and annotated using a search of the SMART (Simple Modular Architecture Research Tool) database (https:// smart. embl. de) [67]. A BLAST search against the NCBI and UniProt databases was performed to find homologous proteins from the common bean G19833 reference genome (Phaseolus vulgaris), soybean (Glycine max), wild soybean (Glycine soja), lupine (Lupinus albus), wheat (Triticum aestivum), maize (Zea mays), rice (Oryza glumipatula), barley (Hordeum vulgare subsp. vulgare), Medicago truncatula, Arabidopsis thaliana and others. Phylogenetic analysis of the candidate genes from common bean was performed by comparison with homologous sequences from other crops. Phylogenetic trees were constructed using the neighbourjoining method as implemented in MEGA 6.0 software with 1000 bootstrap replicates [85]. MEME software [86] was used to identify the motifs of conserved sequences.

Candidate gene expression analysis by qRT-PCR
We selected dry seeds and fresh leaves of BKB and LYD3 for quantitative real-time PCR (qRT-PCR) detection. Total RNA was extracted from the seeds using an improved method combining the cold phenol method [87] and trizol reagent [88], and reverse transcribed into cDNA using a reverse transcription kit (Tiangen, Beijing, China). qPCR was performed in a 20 μL volume consisting of 2 μL of cDNA (1 ng), 0.4 μL of each gene-specific forward and reverse primer (10 μM) (Table S5), 10 μL of 2 × TransStart ® Top Green qPCR Super Mix (+dyeII) and a suitable amount of nuclease-free water. The reaction was run on an ABI 7500 real-time PCR machine (Bio-Rad Corporation, USA) using the following programme: 94 °C for 30 s, 40 cycles of 94 °C for 5 s and 56 °C for 15 s and 72 °C for 10 s. The Actin gene of common bean was used as a reference gene [89]. The 2 −ΔΔCT method [90] was used to analyse the gene expression data from three biological replicates and three technical replicates.