ddRAD sequencing-based identification of inter-genepool SNPs and association analysis in Brassica juncea

Abstract

Background

Narrow genetic base, complex allo-tetraploid genome and presence of repetitive elements have led the discovery of single nucleotide polymorphisms (SNPs) in Brassica juncea (AABB; 2n = 4x = 36) at a slower pace. Double digest RAD (ddRAD) - a genome complexity reduction technique followed by NGS was used to generate a total of 23 million paired-end reads from three genotypes each of Indian (Pusa Tarak, RSPR-01 and Urvashi) and Exotic (Donskaja IV, Zem 1 and EC287711) genepools.

Results

Sequence data analysis led to the identification of 10,399 SNPs in six genotypes at a read depth of 10x coverage among the genotypes of two genepools. A total of 44 hyper-variable regions (nucleotide variation hotspots) were also found in the genome, of which 93% were found to be a part of coding genes/regions. The functionality of the identified SNPs was estimated by genotyping a subset of SNPs on MassARRAY® platform among a diverse set of B. juncea genotypes. SNP genotyping-based genetic diversity and population studies placed the genotypes into two distinct clusters based mostly on the place of origin. The genotypes were also characterized for six morphological traits, analysis of which revealed a significant difference in the mean values between Indian and Exotic genepools for six traits. The association analysis for six traits identified a total of 45 significant marker-trait associations on 11 chromosomes of A- and B- group of progenitor genomes.

Conclusions

Despite narrow diversity, the ddRAD sequencing was able to identify large number of nucleotide polymorphisms between the two genepools. Association analysis led to the identification of common SNPs/genomic regions associated between flowering and maturity traits, thereby underscoring the possible role of common chromosomal regions-harboring genes controlling flowering and maturity in Brassica juncea.

Background

Brassica juncea commonly known as Indian mustard is an important oilseed crop in Indian subcontinent, northern China and eastern European countries. It is widely and extensively grown for seeds which yield an essential oil and condiment; however its young leaves are also used as vegetables or mixed with other salad greens. Brassica juncea has two diverse genepools: the Indian and the east European genepool (exotic) [1]. The east European genepool shows more diversity at the molecular level and has more yield potential while the Indian genepool has narrow genetic diversity with low yield potential [2, 3]. In spite of the two morphological diverse pool, the crop experienced narrow genetic base that might be due to complex allotetraploid genome and domestication [4]. This narrow genetic base has hindered the process of germplasm enhancement as it reduces the chances of finding the diverse alleles of important agronomic traits for their introgression into elite germplasm [5].

The genetic enhancement can be achieved by the transfer of alleles between exotic (European) and Indian genepools using either traditional plant breeding approaches or marker-assisted selection (MAS). While MAS require the identification and use of closely and tightly linked molecular markers with the trait of interest, association mapping does not need prior molecular mapping information and serves as an important tool to identify marker-trait associations on the basis of linkage disequilibrium (LD) only. Association analysis infers significant marker-trait associations by accounting for co-segregated (or co-transmission) alleles at different locations in a genome across a diverse set of mapping population [6], allows fine mapping of traits when used with a dense set of molecular markers. In oilseed Brassica spp. (B. juncea and B. napus) different types of molecular markers were employed with a combination of various models (GLM, Q, PCA and K) to figure out close relationship between various traits and markers. In most of the association mapping studies, the SSR markers were used for population studies due to their usefulness in population genetics inferences and these being highly informative when compared to biallelic markers [7]. However, the high heritability of SNPs also makes them an excellent indicator of genetic diversity and phylogeny in crop species with ancient genome duplications, such as in B. juncea. Various SSR-based genome-wide association mapping studies were conducted in B. juncea [8, 9] and B. napus [10,11,12] for various agronomically important traits. Moreover, Single Nucleotide Polymorphisms (SNPs) are also preferred for fine mapping studies, as reported in B. napus [13,14,15,16,17]. However no such study has been reported in B. juncea, mainly due to non-availability of SNP markers.

The discovery of SNPs in B. juncea has proceeded with a slower pace mainly due to its narrow genetic base, complex allotetraploid genome and highly repetitive regions [18, 19]. The presence of two sub-genomes (A and B) makes SNP discovery and genotyping more difficult and troublesome due to the presence of both homologous and homoeologous DNA sequences. The process of SNP discovery is further complicated with duplications and triplications of A and B genomes due to polyploidization events [20]. To reduce the complexity of genomes, various genome reduction methods are available that uses a set of restriction enzymes and a particular selection process to sequence only the selected set of restriction fragments from multiple genotypes so as to do both the SNP discovery and genotyping at the same time. Advances in the bioinformatics software also support the rapid identification of true SNPs in the individuals.

In this study, a modified ddRADseq approach was followed to partially sequence genomes of six genotypes (three each from Indian and European genepool) of B. juncea for SNPs identification and genotyping. A bioinformatics pipeline was developed using tools available within CLC Genomics workbench for the detection of SNPs (Additional file 1: Figure S1). These SNPs were then used to assess the levels of molecular diversity and population structure among diverse set, and association mapping to identify significant marker-trait associations for six morphological traits.

Results

ddRAD-library preparation and sequencing

The microfluidics-based electrophoresis analyses of pooled library revealed that majority of fragments were represented in the range of 300 to 400 bases following size selection during library preparation (Additional file 1: Figure S2). The sequence-based barcoding followed by pooling and sequencing of six genotypes on Illumina HiSeq 2000 platform generated a total of about 23 million paired-end reads with an average of 3.83 million reads per genotypes. The mean read quality (Phred score) of six samples was 35.02 and about 89% reads had a Q score > 30, indicating that most of the raw data were of good quality. The mean quality score of read 1 (R1) was slightly better than mean quality score of read 2 (R2) (Additional file 1: Table-S3). A slightly better quality of read 1 than 2 was attributed to the fact that the clones within each cluster in a flow cell had least damage due to repeated flushing of flow cell. This difference in quality scores between two reads is also attributed to phasing errors [21].

Sequence analysis and SNP identification

After trimming for low quality sequences, the processed reads were assembled into contigs followed by their alignment to the reference genome of B. juncea (GenBank: LFQT00000000) using default parameters and about 92–94% reads were mapped to the reference genome. The mapping percentages for individual samples ranged from 92.56 (Zem 1) to 94.35 (Donskaja IV). Out of all the mapped reads, more than 80% reads mapped uniquely to a single locus (Table 1).

Table 1 Summary of SNPs obtained in different genotypes

The alignments of contig sequences to the reference genome were used for the identification of SNPs using Probabilistic variant detection method. After filtering for homoeologs, a total of 10,399 single nucleotide variants with a depth of at least 10 reads were found to be distributed among six genotypes.

SNPs in hyper-variable regions and protein prediction

The stringent condition followed during size selection and SNPs identification has led to the retrieval of less number of SNP markers in the genotypes. Although the SNPs were distributed on all the chromosomes there were regions in various chromosomes with high frequency of SNPs as compared to the other regions referred to as hyper-variable regions or SNP hotspots. In all, a total of 44 hypervariable regions or hotspots of SNPs were found on all chromosomes except A08 and B03. Total number of hypervariable regions on these chromosomes ranged from one to five and total number of SNPs in these regions ranged from four to eleven. A BlastX analysis of the hyper-variable sequences identified that nearly 93.2% of these hypervariable regions found to be part of coding sequences (Table 2).

Table 2 Hypervariable regions on different chromosome of Brassica juncea

Morphological analysis of diverse genepools

The diverse core set of B. juncea consisting of 80 genotypes was characterized for various growth and yield traits under two locations in 2015–16 and 2016–17. The data collected over two locations was used to calculate mean values for individual genotypes for days to flowering (DTF), days to maturity (DTM), plant height (PH), siliqua length (SL), seeds per siliqua (SPS) and thousand seed weight (TSW). An analysis of chart correlation for various traits indicated that except DTF and DTM, all the traits were normally distributed in the diverse core set (Fig. 1). In view of the bi-modal distribution for DTF and DTM, the average values for these two traits were used to classify the diverse core set into two genepools namely European (Exotic) and Indian genepools. The individual chart correlation for two genepools indicated normal distribution for all the traits including DTF and DTM.

Fig. 1

Correlation chart for six morphological traits: (a) using 80 genotypes of Brassica juncea of Indian and European genepools indicating bi-modal distribution for DTF and DTM (b) of genotypes of Indian genepool only, and (c) of genotypes of European genepool only

Among two genepools, the traits had divergent correlations with other traits. In Indian genepool, DTF had high correlation with all traits except SPS; while DTF in European genes had significant correlation with DTM only. The DTM in Indian genepool was significantly correlated with PH (0.43) and SL (− 0.43); however the same trait did not show any significant correlation with PH and a positive significant correlation with SL (0.36) and SPS (0.31). PH was significantly correlated with DTF and DTM in Indian genepool but not in European. SL was negatively correlated with both DTF and interestingly with DTM; and positively with TSW in Indian genepool. However, SL had significant but opposite correlation with DTM in European genepool. SPS was not correlated with any of the yield traits in Indian genepool, but was significantly correlated with SL and DTM. In both the genepools, TSW was not significantly correlated with any of the traits but SL in European. The t-test for means for two genepools indicated that the average values for two genepools were significantly different. The p-values for Student’s t-test indicated that the difference in mean values of all the traits among two genepools was highly significant (Table 3).

Table 3 Average values of important traits and their p-values for significance of difference of means between two genepools

Diversity analysis and population structure using SNP markers

A total of 61 SNP markers widely distributed across the B. juncea genome were used for the characterization of core set to develop diversity profile of 80 genotypes. Out of 61 markers, 48 SNP markers were found to be polymorphic. Due to biallelic nature of the marker, a total of 98 alleles were amplified (Table 4). The minor allele frequency ranged from 0.00 to 0.46 with an average of 0.16. The gene diversity and heterozygosity also identified a remarkable degree of variability among the genotypes. The gene diversity value ranged from 0.013 to 0.49 and heterozygosity value ranged from 0.012 to 0.69 with an average of 0.16. PIC (Polymorphism Information Content) values in the present study were found to have ranged from 0.012 to 0.371 with an average of 0.19.

Table 4 Summary of SNP markers used for genetic diversity analysis

The population structure of 80 genotypes was estimated under the Hardy-Weinberg Equilibrium by using STRUCTURE V2.3.4 software. Based on the maximum likelihood and delta K (ΔK) values, the number of optimal groups was identified as two (Fig. 2). A dendrogram constructed using marker allelic data also grouped 80 genotypes into two distinct clusters and the local selection from Turkey forms a separate group. All 80 genotypes were grouped into three major clusters in which cluster I, II and III each contained 29, 50 and 1 genotypes, respectively. Cluster I and II also shows the grouping of genotypes into sub-clusters (Fig. 3). The clustering indicated the ability of SNP markers to group together the related genotypes from a geographical region with high level of accuracy. Cluster I consists of genotype mostly form Indian subcontinent and cluster II consists of exotic genotypes. However, some of the exotic genotypes (EC287711, EC206712, EC491584, EC699038-II and EC699059) were grouped along with the Indian genotypes which may be due to the fact that the allelic composition among these genotypes was identical at some of the loci that were considered in the present study. It may be possible to further refine their grouping patterns by characterization them at greater number of genetic loci as compared to the small subset of 61 SNP markers.

Fig. 2

STRUCTURE analysis indicated genotypes grouping into two sub-populations based on membership coefficients indicated on vertical coordinate

Fig. 3

UPGMA dendrogram showing genetic relationship among 80 different Brassica juncea genotypes. Blue labeled genotypes are of Indian subcontinent, red labeled genotypes are of exotic origin and green labeled genotype is from Turkey

Association mapping analysis

The association analysis to identify markers associated with six yield contributing traits was carried out using a set of 61 SNPs uniformly distributed across the all the chromosomes of B. juncea. In order to determine the true marker-trait associations, we used both p values and marker r² value for association and only those significant associations are considered where the p values were < 10^− 6. Out of 61 SNPs used, 18 SNPs were involved in 45 significant marker-trait associations for six different traits (Table 5). These associations were localized on 11 out of 18 chromosomes of B. juncea with a total of 23 marker-trait associations of A-genome and 22 of B-genome chromosomes. A highest of 16 associations were found for DTF, followed by 13 for SL, 12 for DTM, 2 for TSW and 1 each for SPS and PH. Almost all SNPs, except two – one each on A07 and B02, were involved in multiple associations with different traits. A lone SNP marker on B04 was found to be associated with four different traits; eight SNPs were found to be associated with three different traits followed by associations of six SNPs with two traits each. The SNP markers involved in associations with DTF, DTM and SL, were distributed on both A- and B-genome chromosomes, and these associations were found on multiple chromosomes. The SNPs for SPS, PH and TSW were found to be distributed on single chromosomes only of B-genome. The p-value for all the associations was less than the threshold value as determined by p-value (0.015) of false discovery rate. The p-values for all the associations ranged from 1.26E-05 to 1.15E-18 and the phenotypic variance contribution (r²) ranged from 0.20 to 0.89.

Table 5 Summary of marker-trait associations using SNP markers

Discussion

A plethora of molecular marker-based studies have led to a greater understanding of the genetic make-up of Brassica species. SNP markers have been vital for the (fine) mapping of genes of agronomic importance with the goal of implementing marker-assisted breeding of elite crop cultivars. SNPs are distributed far more frequently in a genome and have been used to develop high-density molecular genetic maps and fine mapping of a region of interest. The abundance of SNPs in genome, low mutation rate and high heritability offsets the disadvantage of bi-allelism. SNPs are found randomly distributed throughout the genome in both repetitive and non-repetitive regions, however those present in the genic/non-repetitive regions are of keen importance. The presence of orthologous regions among the progenitors of allopolyploid genome adds an extra layer of genome complexity in addition to repetitive elements. However, recent advances in reducing the genome complexity coupled with NGS technologies have been highly successful to develop genome-wide SNPs in crops.

In the current study, a pair of restriction enzyme digestion (MseI and SacI) was used for ddRAD sequencing of unique regions of B. juncea. The similar technique of genome complexity reduction has also been employed in several crops [22,23,24] animal [25] and insects [26, 27] species. A number of modifications of this technique have been proposed. In case of other polyploid crop (cotton), GR-RSC (Genome Reduction-Restriction Site Conservation) technique was followed and a combination of EcoRI and BfaI restriction enzymes were used with a size selection between 450 and 600 bp [28] while another study preferred to use a combination of EcoRI and MspI with size selection around 200–400 bp [25].

Following sequencing of genotypes, a total of 2300 MB paired-end sequence data were obtained from six B. juncea genotypes with an average of 383.33 MB from each genotype. Similarly an average of 147.3 MB data was obtained following dd-RAD sequencing of rice [29]. Considering the genome size of B. juncea of 955 MB and the single read sequencing data from six genotypes of 1150 MB, the individual genotype represent an average of 20% of the whole genome and thus, reducing the genome complexity by nearly five folds. Another study on Brassica species reported a reduction of nearly similar genome portion following ddRAD [30]. The mean quality score for both reads ranged from 34.63 to 35.40 and 90% sequence data with a Q score of at least 30 indicated that the sequencing reads were of high quality for reference genome alignment and SNP identification. Similar quality scores for high throughput sequencing runs have been reported with different genome complexity reduction method (SLAF-seq) in tea [31]. Due to high Q score, a large proportion (nearly 83%) of sequence reads were mapped to unique positions in the reference genome indicating the utility of ddRAD method to target unique regions in a genome. The mapping of reads to unique regions also ensured that the SNPs from duplicated or paralogous regions are excluded for further analysis.

Typically, the SNPs are distributed throughout a genome and the average frequency of distribution of SNPs has been found to be between 100 nt to 500 nt. In the present study, the occurrence of 93% of hypervariable regions (hotspots) of SNPs in the coding regions of Brassica juncea with SNPs distributed in upstream, downstream and in the intergenic regions of the coding regions. Most of these hypervariable regions had SNP frequency of less than 10 nt. Further, the detection of 40 genes/coding sequences in the chromosomal regions harbouring SNP hotspots might point to a possible regulatory role of these SNPs in the expression of these genes. Although, few previous studies have reported such SNP hotspots in repetitive regions mostly due to errors of DNA polymerase resulting in strand slippage and unequal exchange [32, 33] or due to presence of mutational hotspots or recombination hotspots [34]. The SNP hotspots along each chromosome were found to be distributed randomly and the number of SNPs involved in such hotspots ranged from four to eleven within 50 nt of chromosomal region in the current study. The role of high selection pressure due to environmental stress could lead to the accumulation of mutated allelic sites in the genic regions that improve survival of the crop under adverse environmental conditions [35, 36].

The high proportion (97%) of functional SNPs across a set of highly diverse genotypes indicated the accuracy of ddRAD technology to invariably target same locus across different individuals during the library preparation and partly due to the improved bioinformatics tools for sequence mapping and SNPs identification for complex and polyploidy crops. The SNPs identified through RAD-seq and its modifications in the previous studies have shown similar functionality levels in other crops as well [37, 38]. The biallelic data obtained from a subset of 61 functional SNPs in the present study was able to group diverse B. juncea genotypes into two major clusters- Indian and Exotic (European) genepool. The diversity and clustering results are in agreement with the previous studies based on SSR and other marker system. The SNP-based diversity analysis also concluded that a small subset of uniformly distributed SNPs would be highly useful for various genetic analyses.

The morphological characterization of six traits revealed very interesting patterns on correlation matrix. The bimodal distribution for DTF and DTM upon combined analysis of all the genotypes indicated that these two traits are controlled by different set of genes in Indian and European genepools. The European genepool has traditionally been domesticated under low-temperature short-day conditions while the Indian genepool is more conducive for sowing in moderate to low temperature conditions found mostly in the north-western plains of Indian subcontinent. The hypothesis of different set of genes controlling DTF and DTM in Indian and European genepools got further strengthened upon getting a unimodal distribution for DTF and DTM in correlation matrices individually for Indian and European genepools. However, the detailed interaction between the genotype and phenotype could be studied by undertaking QTL analysis and other genetic analyses.

In the present study, a common subset of 61 SNPs was used for diversity, population structure and association analyses. For diversity and population structure analyses, the subset of SNPs was able to group 80 genotypes into two distinct clusters, each over-represented by genotypes either from Indian and European (exotic) genepools; which indicated the usefulness of strategy involving the usage of sparse but uniformly localised SNPs for various genetic analyses.

A subset of SNPs representing all chromosomal regions of B. juncea was used to identify significant marker-trait associations. The association analysis using SNP subset was able to localize genes for various agro-morphological trait on different chromosomes, identifying genome regions for undertaking fine mapping of traits/genomic regions with large number of molecular markers. A majority of SNPs identified associations with multiple traits thus essentially indicating either the clustering of genes for multiple traits or involvement of same set of genes regulating multiple traits in the same genomic regions. Among these traits, DTF and DTM had invariably common SNP/genomic region associated with them, thus implying that the genes for these two traits are clustered together and/or likely have correlated/coordinated expression of genes. A recent study, using F₂ mapping population, in Brassica napus has also identified the co-localisation of QTLs (and eQTLs) for flowering time and various growth-related morphological traits to a common genomic region of chromosome A10 [39]. In another study, QTLs for various quality and nutritional traits were again mapped to common regions of a genetic map of a DH (double haploids) mapping population in Brassica napus [40]. Such clusters of QTLs for multiple traits were also reported using chromosome segment substitution lines (cssls) in Brassica rapa [41]. High correlation between DTF and DTM traits, in the current study, among both Indian and European genepools also indicate the high probability of association of common genomic regions (and SNPs) for both the traits as reported in one of the earlier study as well [42].

The presence of a common ancestral genome between three polyploidy species led to the identification and comparison of association analysis results. In the current study, the associations for DTF were mapped to A- and B- genome chromosomes. Similarly genes for flowering time have been identified on both A- and B- subgenomes of B. juncea [43]. Two highly significant associations for DTF were identified each at 6.8 MB (A06_6796237) and 23.4 MB (A06_23478761) in the current study are in agreement with the results for flowering time related (FTR) genes. Thirty three flowering time related (FTR) genes were identified on chromosome A06 between 7.2 MB – 21.6 MB regions using transcriptome analysis [44]. The association analysis results of the current study indicated that a subset of sparse but uniformly localised SNPs would be highly useful to demarcate genomic regions for traits of interest.

Conclusion

This is the first report of use of ddRAD-seq for the development of SNPs in Brassica juncea. The SNPs were developed initially from sequence comparison of six genotypes only; however the SNPs were found to be functional when tested on a diverse set of genotypes. The SNPs used for association analysis were also found to be significantly associated with six morphological traits. Given the fact that Brassica juncea has narrow genetic base, the SNPs identified in the current study would form an excellent source for various genetic studies including linkage mapping, fine mapping and association analysis.

Methods

Plant material and DNA extraction

A set of six B. juncea genotypes (three each from Indian and Exotic germplasm) were selected for use in ddRAD library preparation. Pusa Tarak (BJI-1), Urvashi (BJI-2) and RSPR-01 (BJI-3) were selected from Indian genepool and Zem 1 (BJE-1), Donskaja IV (BJE-2) and EC287711 (BJE-3) were selected from European (exotic) genepool. Seeds were procured from (Dr. Deepak Pental) University of Delhi (South Campus), India and National Bureau of Plant Genetic Resources, New Delhi, India. SNP genotyping was performed on 80 diverse B. juncea genotypes that were procured from Plant Gene Resources, Agriculture and Agri-Food, Canada and Genetics & Plant Breeding Department, SKUAST-Jammu, India (Additional file 1:Table-S4). Total genomic DNA was isolated using modified SGS buffer method [45] and purified DNA was used for dd-RAD library preparation.

Morphological data evaluation and statistical analysis

The phenotypic data of diverse core set of B. juncea was also recorded from two different locations in 2015–16 and 2016–17. The data were collected for six traits: days to flowering (DTF- number of days from sowing to the date when 50% of the plants had their flower opened in each plot), days to maturity (DTM- number of days from sowing to the date when pods on 75% of the plants in each plot were turned browned), plant height (PH- in meters), siliqua length (SL- in centimeters), seeds per siliqua (SPS- average number of the seeds present in single pod/siliqua) and thousand seed weight (TSW- weight in grams of the 1000 seeds collected in random). The traits value of each genotype was defined as an average of two replicates in the same location. The correlation coefficients between traits were determined using Student’s t-test and the variance components were also calculated.

ddRAD library development and NGS sequencing

The ddRAD-seq protocol [25, 46] was used with slight modification for the construction of sequence-barcoded reduced representation libraries (RRLs) from six Brassica juncea genotypes. For ddRAD library preparation, ten microgram of purified DNA was digested to completion with Mse I and Sac I. The digested DNA was separated on 0.8% agarose gel; fragments between 300 and 400 bp were gel excised and eluted. The eluted and purified DNA was then end repaired, short dA-tail was attached and ligated with the adapters following manufacturer protocol. The ligated DNA was amplified using PCR to enrich and add the Illumina specific index and flow cell annealing sequences to the fragmented DNA. For each six genotypes, six different index sequences were used so at to facilitate the process of pooling. All six DNA samples were normalized to a final concentration of 50 ng/μl and pooled to reach a final volume of 300 μl to generate a reduced representation library. The pooled dd-RAD library was then sequenced using Illumina HiSeq 2000 to generate 100 bp paired-end reads.

Sequence preprocessing and SNP detection

The ddRAD-seq reads obtained after sequencing were bioinformatically analyzed using CLC Genomics Software in order to obtain a high quality SNP set. The paired end sequencing reads were subjected to a series of steps (demultiplexing, trimming, mapping with reference genome, local realignment, SNPs detection and annotation with flanking sequences) through a pipeline. The following filtering scheme (Fig. 4) was used to maximize the retention of true genic polymorphic SNPs: (1) trimming of 13 bases from forward and 3 bases from reverse end, (2) mapping parameters were set to- mismatch cost: 2, insertion cost: 3, deletion cost: 3, length fraction: 0.5, similarity fraction: 0.95 and we have selected to perform local alignment instead of global alignment as it allows the ends to be left unaligned if there are many differences from the reference at the ends, (3) probabilistic SNP detection method was used for SNP detection from mapped reads with parameters- minimum coverage: 4, variant probability: 98.00 and ploidy: 2 and (4) flanking sequence of 400 bp. For mapping of reads, Brassica juncea genome was used as a reference genome [47].

Fig. 4

Workflow of various steps involved in SNPs identification

Validation of SNPs and genotyping

A subset of 61 SNP loci was selected with 3–4 SNPs from each chromosome and was validated across the diverse set of B. juncea. The sequences flanking each SNP were used to synthesize forward, reverse and iPLEX universal extension primer using Agena CXassay design suite V2.0 software. The forward and reverse PCR primers were diluted to the concentration of 100 μM, while iPLEX universal extension primers were diluted to the concentration of 500 μM. The experimental procedure included- (1) multiplex PCR using forward and reverse primer, (2) SAP (Shrimp Alkaline Phosphatase) clean up reaction, (3) iPLEX extended reaction with the amplified product, (4) resin cleanup reaction to remove salts, (5) spotting of primer extended product on spectro-chip and (6) spectro-chip detection using MALDI-TOF mass spectrometry. The genotype calls were evaluated through MassARRAY TYPER 4.0 software.

Population and diversity analysis

The SNP genotyping data were used for population structure and genetic diversity analyses [48] following Singh et al. [49]. The posterior probabilities (qK) were estimated with 10,000 burn-ins followed by 100,000 iterations. For structure analysis, the diverse population was assumed to be following an admixture model and correlated allele frequencies with no prior population information. The structure analysis was performed with 5 replicates for each K ranging from 1 to 5. The ΔK was calculated using Structure Harvester software [50] to obtain an optimal value of K. The membership coefficient with a threshold of 70% for each replicate of structure analysis was used to generate a Q matrix using the software CLUMPP [51]; followed by plotting of Q matrix using DISTRUCT software [52]. The polymorphic information content (PIC) value and allele frequencies were calculated using Powermarker v3.51 [53]. The unweighted neighbor joining tree method was implemented in Darwin5 software [54] for constructing a phylogenetic tree; and the bootstrap value for this tree was determined by re-sampling loci at 1000 times.

Gene identification and annotation using database

Flanking sequence of SNPs/ hyper-variable regions were compared against the B. juncea database using BLASTX (cutoff E-value of 1E-10) to identify the corresponding sequences in the protein database [55].

Association analysis

Association analysis was performed by using the genotypic (SNPs) and phenotypic data of the diverse Brassica genotypes and population structure data (Q matrix) by using TASSEL software [56]. Marker–trait association analysis was conducted using TASSEL 3.0 software along with the GLM procedure keeping significant threshold for the association at P < 0.01.