Until recently, precise location of genes and marker assisted selection was long thought in legumes such as blackgram due to lack of dense molecular maps. However, advances in next generation sequencing based high throughput genotyping technologies such as QTL-seq have revolutionized trait mapping in marker-orphan crops. Using QTL-seq approach, we have identified a large effect QTL for resistance to Mungbean yellow mosaic India virus (MYMIV) in blackgram variety Mash114. MYMIV is devastating disease responsible for huge yield losses in blackgram, greengram and other legumes. Mash114, a blackgram variety developed by Punjab Agricultural University, India showed consistent and high level of resistance to MYMIV since last nine years. Whole genome re-sequencing of MYMIV resistant and susceptible bulks derived from RILs of cross KUG253 X Mash114 identified a large effect QTL (qMYMIV6.1.1) spanning 3.4 Mbp on chromosome 6 explaining 70 per cent of total phenotypic variation. This region was further identified as an inter-specific introgression from ricebean. The annotation of the introgressed segment suggested presence of clusters of disease resistance domains of which nine candidate genes were found to have role in virus resistance. Linkage mapping using KASP markers developed from these nine candidate genes identified the 500 kb genomic region equaling 1.9 cM on genetic map linked with MYMIV. The three KASP markers closely associated with MYMIV originated from serine threonine kinase, UBE2D2 and BAK1/ BRI1-ASSOCIATED RECEPTOR KINASE genes. These KASPs can be used for marker assisted transfer of introgressed segment into suitable backgrounds of Vigna species.
Research Article
A large effect QTL introgressed from ricebean imparts resistance to Mungbean yellow mosaic India virus in blackgram (Vigna mungo (L.) Hepper)
https://doi.org/10.21203/rs.3.rs-1759872/v1
This work is licensed under a CC BY 4.0 License
You are reading this latest preprint version
Blackgram
Vigna mungo
ricebean
Vigna
umbellata MYMIV
yellow mosaic disease
QTL-seq
QTL
KASP
candidate genes
Here, we report identification of a large effect QTL conferring Mungbean yellow mosaic India virus resistance introgressed from ricebean in blackgram variety Mash114. The tightly linked KASP markers would assist in marker-assisted-transfer of this region into Vigna species infected by MYMIV.
The genus Vigna and subgenus Ceratotropis are the major source of essential food proteins, predominantly in Southeast Asia. Cowpea (Vigna unguiculata (L.) Walp.), mungbean (V. radiata (L.) R. Wilczek), blackgram (V. mungo (L.) Hepper) and azuki bean (V. angularis (Willd.) Ohwi & Ohashi) are the most important species among the others in this genus. Black gram, also known as urdbean, is a self-pollinating annual crop with diploid chromosome number (2n = 2x = 22) and moderately sized genome of 574 Mbp (Arumuganathan and Earle 1991). Blackgram is nutritionally rich crop as it consists of about 56% carbohydrate, 25% protein, 4% minerals, 2% fat, 0.4% vitamins and high content of phosphoric acid. It is a vital grain legume with digestible proteins, low content of flatulence and also replenishes nutrient status of soil by biological nitrogen fixation (Kumar and Pandey 2020). The main hindrance in unleashing aforesaid benefits of crop is low productivity of blackgram due to biotic (yellow mosaic disease, cercospora leaf spot and powdery mildew) and abiotic stresses (heat, drought and pre- harvest sprouting) among other factors. Amongst the different biotic stresses, yellow mosaic disease (YMD) caused by Mungbean yellow mosaic India virus (MYMIV) is a serious concern causing up to 85–100% yield loss (Bashir et al. 2006). Although, YMD is prevalent in all blackgram growing countries except Australia, but heavy incidence of the disease is mainly reported from India, Pakistan and Bangladesh causing devastating crop losses (Biswas et al. 2008).
MYMIV is an important virus of leguminous crops, transmitted by an insect vector i.e. whitefly Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae). MYMIV belongs to family Geminiviridae, the second largest family of plant viruses responsible for many devastating epidemics (Varma and Malathi 2003) and genus begomoviruses, the largest genera of this family (Qazi et al. 2007). YMD is caused by two main viral species, of which Mungbean yellow mosaic India virus (MYMIV) is reported to be more prevalent in Northern, Eastern India and Central India (Naimuddin and Akram 2010) and Mungbean yellow mosaic virus (MYMV) in Southern India particularly peninsular region (Karthikeyan et al. 2004; Malathi and John 2008). Since the present study was conducted in northern part of country, therefore term ‘MYMIV’ was used for referring the viral pathogen causing YMD. The symptomatology of MYMIV was first studied by Nariani (1960) and its name of yellow mosaic disease is because of appearance of mosaic of yellow and green patches on leaves of infected plants. The most prominent symptom is found on young leaves as small yellow spots on vein-let and ultimately covering the entire lamina. Later on, the spots spread and the complete leaf turns yellow hampering the photosynthetic ability of plants and rendering stunted growth and low yield.
Mapping attempts for MYMIV resistance in blackgram in the past could not be much successful due to lack of availability of molecular markers, dense genetic maps, any cytogenetic stocks and only recent genome sequencing of Vigna mungo (Pootakham et al. 2021). Basak et al. (2005) used resistance gene analog (RGA) markers to map YMD resistance in a resistant mutant of cultivar T9. Souframanien and Gopalakrishna (2006) screened 100 ISSR markers from related species to map YMD resistance in 58 RILs derived from TU 94-2 (resistant) and V.mungo var selvestris. They identified one linked ISSR which was then converted into SCAR marker. Gupta et al. (2013) attempted mapping of YMD resistance from DPU 88 − 31 resistance cultivar by bulked segregant analysis (BSA) using 361 SSR markers from sister crops such as cowpea and mungbean. However, due to lack of high-throughput and dense molecular markers, precise location and molecular markers tightly linked to YMD resistance could not be identified. In addition, most of the resistant sources used in various mapping studies were directly and indirectly derived from same cultivars which suggested narrow genetic base of blackgram.
The advent of next generation sequencing based high-throughput genotyping methods and availability of reference genome sequence have accelerated molecular mapping and map based cloning of genes in most of the economically important crop plants (Tuberosa and Salvi 2007). QTL-seq has emerged as an important and high-throughput approach to identify genomic regions linked with trait of interest (Takagi et al. 2013). It analyses genome wide sequence differences in two contrasting bulks to identify the causal SNPs governing the target phenotype. To quote a few, the approach has been instrumental in identification of major effect QTLs in rice (Takagi et al. 2013), wheat (Wang et al. 2021), tomato (Illa-Berenguer et al. 2015) and squash (Ramos et al. 2020).
For developing durable resistance to MYMIV in blackgram, there is a need to identify novel sources of resistance to MYMIV. Punjab Agricultural University (PAU), India has developed a MYMIV resistant line Mash114 from an inter-specific cross of blackgram and ricebean where MYMIV resistance from ricebean has been transferred to Mash114 using repeated backcrossing and selection (Singh et al. 2013). Mapping and identification of tightly linked markers to MYMIV from Mash114 will facilitate accelerated transfer of novel resistance to elite blackgram genotypes. In the present study, we report a major effect QTL introgressed from ricebean imparting resistance to MYMIV in Mash114. We have also developed KASP markers that could be used for marker-assisted transfer of this trait. In addition, suggested candidate genes could be explored further for map based cloning and functional characterization of gene(s) responsible for MYMIV resistance.
Plant material and screening against MYMIV
Plant material included MYMIV resistant blackgram line, Mash114, MYMIV susceptible line KUG253 and a set of 160 RILs derived from their cross. Mash114 is a released variety of Punjab Agricultural University (PAU), India which was developed from an inter-specific cross of V. mungo cv. Mash 338 and V. umbellata cv. RBL1 (Singh et al. 2013). KUG253 is high yielding advanced breeding line of V. mungo but shows highly susceptible reaction to MYMIV.
The 160 RILs along with parents were planted in three replications in randomized complete block design at PAU, Ludhiana (hotspot for MYMIV in India). Each entry was planted in 2m row with plant × row spacing of 10 × 30 cm and 6 + 3 (entry + infector) rows in one bed. For proper spread of disease, MYMIV susceptible cultivar SML1082 was used as infector row and was planted all around the borders and middle of each bed.
MYMIV disease scoring (DS) was performed on scale of 1-9 (Singh and Bhan 1999) with DS=1 (No visible symptoms), DS = 3 (minute yellow specks covering 1-15 per cent leaf area), DS = 5 (Yellow mottling and discoloration of 15.1-30 per cent leaf area), DS = 7 (Severe yellow discoloration of leaves covering 30.1-75 per cent of foliage, stunting of plants and reduction in pod size), DS = 9 (Severe yellow discoloration of leaves covering 75.1-100 per cent of foliage, stunting of plants and no seed formation). Data was recorded from 20 plants of each progeny and averaged. Disease was scored after the appearance of symptoms on susceptible check and continued at an interval of every 7 days, till the susceptible check became fully susceptible (DS=9).
DNA extraction, bulks generation and whole genome re-sequencing
Genomic DNA was isolated from young leaves of RILs and parental genotypes using CTAB method (Rogers and Bendich 1989). The quality and quantity of genomic DNA was analyzed using Nanodrop ND-2000 spectrophotometer according to manufacturer’s manual (Thermo Fisher Scientific, USA). Based on phenotyping, 40 RILs with DS of 1 and 40 RILs with DS =7 were selected. DNA was extracted from each RIL and was normalized to 1microgram/microlitre concentration. Equal quantity of DNA of each individual RIL was then pooled to form MYMIV resistant bulk (MYMIV-R), and similarly MYMIV susceptible bulk (MYMIV-S). Whole genome of both parental lines and two bulks i.e. MYMIV-R and MYMIV-S were sequenced at 10X using pair-end sequencing on Illumina NovaSeq 6000 sequencer which was outsourced to NGB Diagnostics, New Delhi, India.
Variant identification
The short reads obtained from paired end sequencing of parents and bulks were filtered for low quality bases using FastQC (Babraham Institute, England). The adapter sequences in the reads were trimmed by Trimmomatic (Bolger et al. 2014). Thereafter, the high quality clean reads were aligned to reference genome of blackgram (Pootakham et al. 2021); https://www.ncbi.nlm.nih.gov/data-hub/assembly/GCA_013427195.1/) using BWA (Li and Durbin 2009) with default parameters. The resulting file was saved in SAM (Sequence Alignment Map) format, which was then converted to BAM (Binary Alignment Map) format with the help of Samtools (Li et al. 2009). Further, picard tool (https://broadinstitute.github.io/picard) was used to modify read groups to make them compatible for variant calling followed by sorting the BAM files using Samtools. Variant calling was performed using Genome Analysis Toolkit (GATK, McKenna et al. 2010) to obtain single nucleotide polymorphism (SNP) and insertion and deletions (InDels). The variant file was filtered based on missing data ≥ 0.9 and minimum allele frequency ≥ 0.05 using VCFtools (Danecek et al. 2011). The filtered variants including SNPS and InDels were converted to table file using VariantsToTable format.
QTL-seq analysis
An R package QTLseqr (Mansfeld and Grumet 2018) was used to identify genomic regions linked to MYMIV in two contrasting bulks using read depth of 10, minimum GQ of 30 and refAlleleFreq of 0.20. The ratio of reference allele count (high bulk reference allele depth + low bulk reference allele depth) to total read depth in combined bulks was used to calculate reference allele frequency. Genome wide G statistics (G' value) and SNP index for each SNP site were then calculated from observed and expected allele depth of the reference and alternate allele. The tricube smoothed G prime (G') and delta SNP-index {(Δ) SNPindex} values of each SNP were estimated using two main functions viz., runGprimeAnalysis and runQTLseqAnalysis respectively (Magwene et al. 2011). The p-values were estimated for each SNP using non-parametric estimation of the null distribution of G’ statistics. Tricube-smoothed (Δ) SNP-index was calculated using a 1Mb sliding window. The analysis was bootstrapped 10000 times and two sided 95 and 99 per cent confidence intervals were calculated for each Δ (SNP-index) using script runQTLseqAnalysis. The plotQTLStats function was used to plot graphs for number of SNPs, (Δ) SNP-index, G′ value and -log10(p-value). Finally, getSigRegions and getQTLTable functions exported the regions considered significant at the specified confidence interval to identify putative QTL(s).
Annotation of candidate region
The reference genome annotation of blackgram was not available; therefore Augustus 3.2.1 was used to predict coding sequences (CDS) and gene prediction in the QTL region by ab initio gene predictor (Stanke et al. 2006). CDS were aligned to NCBI nr database using BLASTx with eukaryotic genes and annotated functionally using OmicsBox 1.3.11 (https://www.biobam.com/downl oad-omicsbox/) (evalue ≤1e-5). CDS were subsequently mapped to retrieve Gene Ontology (GO) terms and further mapping of enzyme code with KEGG (Kyoto Encyclopedia of Genes and Genomes) database (http://www.genome.jp/kegg). In order to enrich the knowledge of gene functions and their precise locations in blackgram genome, well saturated genome annotation of a related crop species (V. unguiculata) was used from NCBI (Muñoz-Amatriaín et al. 2017).
KASP development and linkage mapping
KASP (Kompetitive Allele Specific PCR) were designed from SNPs in introns and exons of candidate genes spanning the target region using Primer 3.0 (Additional file 1: Table S1 and Table S2). KASP assays were performed in a 384-well by following instructions of LGC Genomics (LGC, Middlesex, UK) with PCR profile: 94 °C for 15min, 95 °C for 20 s, 55 °C for 60s (−1 °C/cycle, 10 cycles in total) and 94 °C for 10 s; 57°C for 60 s (30cycles). The reaction mixture and concentration of DNA was prepared as per manufacturer’s protocol. The Synergy H1 fullfunction microplate reader (FLUO star Omega, BMG Labtech, Germany) was utilized to read the plates and generate fluorescence signal. Linkage groups were constructed using ICIMapping 4.2 software, by grouping module at LOD >3.0 to interpret linkage of KASP markers. Kosambi map function was used to compute genetic distances, single marker analysis and inclusive composite interval mapping (ICIM) was performed by ICIMapping 4.2 (http://www.isbreeding.net) with threshold LOD ≥ 3.0 and 0.05 as significance level (Churchill and Doerge 1994) to identify QTL location.
Screening of parents and RILs against MYMIV
The resistant (Mash114) and susceptible (KUG253) parents showed clear phenotypic differences (Figure 1) while F1 showed resistant reaction (DS=1) indicating dominance of resistance over susceptibility. The frequency distribution in RILs showed symmetric distribution (Figure 2). There were 70 asymptomatic RILs, thus scored as resistant (mean DS=1) while 90 RILs were MYMIV susceptible showing characteristic yellow mosaic symptoms. This distribution of resistant (R) and susceptible (S) fitted the RIL population in 1R:1S ratio (ꭓ2=3.5, p ≥ 0.05) indicating single dominant gene conferring MYMIV resistance. Of the 90 susceptible RILs, 50 RILs showed very high chlorosis in leaves with stunted growth thus given the score of DS=7, whereas 17 RILs showed low chlorosis and scored as DS=3, while 23 RILs were scored DS=5 with moderate chlorosis in leaves and retarded growth, when compared with MYMIV score of infector row (mean DS=7) (Figure 2). Based on MYMIV reaction data, 40 RILs showing extreme resistance and 40 RILs showing extreme susceptible reaction were selected to constitute MYMIV-R and MYMIV-S bulks respectively.
Whole genome re-sequencing analysis and variant detection
A total of 37.6 Gbp clean reads with QC > 30 were obtained from re-sequencing of two bulks (MYMIV-R, MYMIV-S) and two parental genotypes (Table 1) that represented 90 per cent of the original sequencing data. An average sequencing depth of 16.5 x with aggregate 212.9 million clean reads were obtained after initial quality check and mapping of these clean reads to the CN80 blackgram reference genome gave a mean mapping coverage of 92 per cent. GATK identified ~18 million genetic variants encompassing 1522432 SNPs and 301425 InDels in susceptible parent ‘KUG253’ whereas, 1674446 SNPs and 335549 InDels were identified in the resistant parent, Mash114, with a total of 2795896 variants. Among 2795896 variants, transitions (2017687) outnumbered transversions (666488) for both the parental genotypes. All InDels were filtered out for simplicity of analysis and retained after final analysis. SNPs were filtered with maximum missing 0.9, this led to a total of high quality 1861060 SNPs available for QTLseqr analysis, covering 11 chromosomes (492 Mbp of total 570 Mbp genome). However, the SNPs present on smaller contigs were excluded for further analysis.
QTL-seq analysis and identification of candidate QTL region
A total of 588821 high quality SNPs after filtering by QTLseqr default parameters (Additional file 1: Figure S1) were used for QTL-seq. To infer the QTL region conferring MYMIV resistance, the genome wide comparison SNP-index of MYMIV-R and MYMIV-S bulk over the entire length of genome was performed using QTLseqr. The value of SNP-index equaling ‘0’ implies that all the sequencing reads at a particular genomic position in a bulk originated exclusively from susceptible KUG253 genome, while SNP index value of ‘1’ implied that reads originated entirely from resistant Mash114 genome. The SNP-index values of MYMIV-R and MYMIV-S bulks identified four QTLs named as qMYMIV6.1, qMYMIV6.2, qMYMIV6.3, qMYMIV6.4 to be associated with MYMIV resistance on chromosome 6 harboring QTL region containing 16613, 2046, 7014 and 723 SNPs, respectively (Table 2, Fig. 3). The physical distances between the QTLs qMYMIV6.1, qMYMIV6.2, qMYMIV6.3 and qMYMIV6.4 were 2, 5 and 4 Mbp, respectively. The QTL region qMYMIV6.1 was the largest spanning over 8.30 Mbp, followed by qMYMIV6.3 with length of 4.93 Mbp, qMYMIV6.2 with 2.0 Mbp and qMYMIV6.4 covered 0.52 Mbp on chromosome 6.
The 3.4 Mbp region of QTL qMYMIV6.1 from 5.6 to 9 Mbp of chromosome 6 was characterized by average SNP-index lower than 0.2 (lowest value = 0) and higher than 0.8 (highest value = 1) and temporarily designated as qMYMIV6.1.1. While, the rest of the genomic regions of qMYMIV6.1 and remaining three QTLs did not show similar trend; therefore not considered for further analysis. In the qMYMIV6.1.1 region, all the component RILs of MYMIV-R bulk had resistant parent Mash114 type allele and all component RILs constituting MYMIV-S bulk had susceptible parent KUG253 type allele as evident from individual reads in the bulks on integrated genome viewer (IGV) (Figure 4). This strongly suggested that qMYMIV6.1.1 could be probable region for gene underlying MYMIV resistance. Moreover, peak ∆(SNP-index) of 0.65 and peak Gprime value of 42.08 at 8.7 Mbp position of qMYMIV6.1.1 further supports this as putative QTL region. Therefore, qMYMIV6.1.1 was selected for further analysis.
Annotation and identification of candidate genes
A total of 1167 CDS were predicted in qMYMIV6.1 region which were annotated by sequence similarity searches using BLASTx algorithm against NCBI non-redundant (nr) protein database of which 46 % of cds had significant similarity to already known proteins while 53 % had no hits in nr database and referred to as uncharacterized protein (Additional file 2: Data S1). The functional annotation by Blast2GO assigned functional terms and 9330 gene ontology numbers to 547 CDS differentiating their molecular function, biological process and cellular component. Of these, 150 cds mapped to KEGG pathways having role in disease resistance, were selected. The disease resistance genes present on targeted 3.4 Mbp QTL region of MYMIV6.1.1 were visualized on IGV. The 20 genes were identified with read depth more than 40 and in silico polymorphism in parental genotypes as well as corresponding bulks. Of these 9 candidate genes were selected based on GO function, enzyme activity and KEGG pathway for potential involvement in virus resistance.
KASP marker development, linkage mapping and identification of major effect locus
To further validate the association of qMYMIV6.1.1 region to MYMIV and understand genetic relationship of SNPs in the region, a total of nine SNPs present in both introns and exons of selected candidate genes and spanning the QTL region, were extracted. The selected SNPs were converted to KASP markers (Additional file 1: Table S1 and Table S2). The nine KASP markers were used to genotype 160 RILs (Fig. 4) that generated a genetic linkage map of 15.98 cM. KASPs were closely placed on genetic linkage map such that the genetic distances between adjoining markers VM608, VM601, VM609, VM604, VM605, VM610, VM602, VM606 and VM607 were 3.73, 2.76, 1.58, 2.57, 0.78, 1.12, 2.06 and 1.38 cM, respectively. All the KASP markers showed linkage to MYMIV resistance, Further, inclusive composite interval mapping mapped a large effect QTL at maximum LOD value of 42.4 showing PVE = 70.89 % and additive effect = -2.3201 associated with MYMIV. This major effect locus spanned a genetic distance of 1.9 cM including three markers viz., VM605, VM610 and VM602. The interval between markers VM605 and VM610 (0.77 cM) showed LOD value of 42.4 at P > 0.05 and the interval between markers VM610 and VM602 (1.12 cM) showed LOD value of 41.9 at P > 0.05 (PVE = 71.28 %, Additive effect = -2.3201). Thus, it was confirmed that VM602, VM605 and VM610 were tightly linked to MYMIV resistance gene spanning a genetic distance of 1.9 cM which corresponds to 0.5 M physical region on V. mungo genome.
MYMIV resistance in Mash114 introgressed from V. umbellata
Genotyping of V. umbellata genotype RBL1 using KASP markers indicated 2.1 Mbp fragment on chromosome 6 from genomic position 6978055 bp (VM609) to 9085008 bp (VM605) in Mash114 had same alleles as RBL1 (V. umbellata) (Fig. 6). These findings intrigued convincing clue that the gene(s) controlling MYMIV resistance in V. mungo cv. Mash114 are present on introgressed segment from V. umbellata. The candidate region (0.5 M) of introgressed segment harboured 80 genes as predicted by Augustus (Additional File 2: Data S2), of which 9 genes were involved in disease resistance pathway carrying 31 non-synonymous SNP mutations.
In addition, comparison of physical and genetic positions of KASP markers in blackgram revealed similar order of arrangement in both maps (Fig. 6) with few exceptions. Since, genome of cowpea (Vigna unguiculata) was well annotated than blackgram (Vigna mungo), sequence similarity of V. unguiculata with V. mungo was inferred to enrich the annotations of V. mungo and fill the gaps. It was found that the QTL region (0.5 Mbp) mapped by KASP markers on chromosome 6 of blackgram was syntenic with 430 Kb genomic fragment of chromosome 2 in Cow pea (Vigna unguiculata). This syntenic fragment (24-26Mbp) of chromosome 2 of V. unguiculata was further disintegrated to reveal the component genes. Interestingly, abundance of plant disease resistance gene clusters including 37 domains of LRR receptor like serine threonine kinases and thaumatin like protein, domains of putative disease resistance protein RGA, disease resistance RPP13-like protein 4, F-box protein , E3 ubiquitin-protein ligase in QTL region.
Our study has identified a large effect QTL associated with MYMIV resistance in V. mungo and further identified tightly associated KASP markers that can be used in breeding program. However further narrowing down of this region will help to identify putative candidate gene(s) responsible for MYMIV resistance and unleashing its function thereof.
Indian farmers are restraining from the cultivation of Vigna mungo alongwith with other Vigna crops, despite its multiple benefits to human & soil health and its shorter duration of about 50–60 days. This is due to low yields in this crop, as not much efforts are invested for developing farmer friendly varieties of Vigna in comparison to emphasis given to major cereal crops and some major pulses. MYMIV transmitted by whitefly, is a serious concern for low yields of crop which should be addressed meticulously (Project Coordinators Report 2018; Mishra et al. 2020). Yield losses from MYMIV may reach up to 85–100% when conditions are suitable for disease development (Bashir et al. 2006). Since, domestication and genetic erosion caused by human intervention with recurrent use of genotypes belonging to common ancestry has led to narrow genetic base in blackgram. Thus, a systematic program for the identification, introduction and mapping of MYMIV resistance genes is required which may help to infuse the diverse and novel genes into high yielding germplasm to broaden the genetic base of crop. Wild germplasm and related species are goldmine for any crop species, and these serve as reservoir of multitude novel alleles/genes of economic importance (Pandiyan et al. 2008). The MVMIV resistance derived from V. umbellata in the blackgram variety Mash114 developed by PAU, is maintaining its resistance since last 9 years indicating the effectiveness of resistance derived from related species (Singh et al. 2013). Mapping of MYMIV resistance from Mash114 and identification of tightly linked markers would facilitate the marker-assisted transfer of this trait in elite genotypes.
QTL-seq is an NGS based high throughput genotyping technique that identifies genomic region associated with target trait which is challenging task in crop like blackgram through traditional gene/QTL mapping limited by unsaturated linkage maps. The approach is a shorter and cost effective route to scan whole of genome for identification of causal variants. Whole genome resequencing of parents identified more genetic variants in resistant parent than susceptible against reference genome as anticipated from inter-specific lineage of Mash114. The comparison of SNP-indices along with G’ and p-value statistics of MYMIV-R and S bulks covering the entire blackgram genome identified four consecutive QTL regions (qMYMIV6.1, qMYMIV6.2, qMYMIV6.3 and qMYMIV6.4) on chromosome 6 controlling MYMIV resistance. The four QTL regions illuminated pattern of recombination over entire length of chromosome 6 as MYMIV linked regions with no recombination were detected as QTLs peaks and other unlinked regions were dissolved by recombination events of RILs as a result of pooling of DNA of R and S RILs used for developing bulks. Among the four QTLs, qMYMIV6.1.1 a sub-fragment of qMYMIV6.1 was selected as a strong candidate region harbouring MYMIV resistance gene for further analysis, because it displayed mean ∆(SNP-index) of 0.25, peak ∆(SNP-index) of 0.65 and peak Gprime value of 42.08 at 8.7 Mbp position. The other three QTL-regions (qMYMIV6.2, qMYMIV6.3 and qMYMIV6.4) were rejected for further analysis as these regions harboured reads derived from both resistant and susceptible parent in both bulks as viewed in alignment file on IGV. In contrast, region spanning qMYMIV6.1.1 in MYMIV R bulk showed reads exclusively derived from R or S parent.
The annotations of blackgram predicted in present study were enriched using annotations from cowpea (V. unguiculata), because it was well-annotated Vigna species. The qMYMIV6.1.1 was found syntenic with a fragment on chromosome 2 of cowpea genome, although the orientation of chromosome was inverted in cow pea. Interestingly, the annotations from blackgram and cow pea genome suggested presence of clusters of disease resistance genes in this region e.g. 37 domains of LRR receptor like serine threonine kinases, putative disease resistance protein RGA, disease resistance RPP13-like protein 4, F-box protein and E3 ubiquitin-protein ligase etc. The predominance of clusters of disease resistance genes in the candidate region signaled profound network of resistance genes embedded on chromosome 6. We selected 9 candidate genes for KASP marker development with defined role in host defiance against viruses. The KASP derived from candidate genes, i.e. VM608 from transcription initiation factor IIB, (GO:0019083), VM601, VM602 (UBE2D2), VM606, VM607, VM610 (BAK1/ BRI1-ASSOCIATED RECEPTOR KINASE) from serine/threonine-protein kinase (GO:0004674), VM609 from IRES-dependent viral translational initiation (GO:0075522), VM604 from Serine/threonine/tyrosine kinase activity (GO:0003779), VM605 from plant NLR complex (GO:0042742) elucidated meiotic cross-overs and linkage pattern with MYMIV resistance gene.
Protein serine/threonine kinase (STKs) class of domains was predominantly present in QTL region qMYMIV6.1.1 and the same was targeted in current study for designing KASPs. STKs is a class of receptor like kinases (RLKs), evolutionary conserved proteins having role in plant defense (Afzal et al. 2008) as reported for powdery mildew in wheat (Cao et al. 2011), stripe disease resistance in rice (Lee and Kim 2015) and northern corn leaf blight resistance (Hurni et al. 2015). Plant proteins conferring resistance against pathogens share similar domain, motif, and structure. Such a large-scale similarity is derived from evolutionary conservation indicates that similar signaling pathways underlie plant defense cascades. The recent reports indicated role of STKs in plant defense against viruses also. Geminivirus’s nuclear shuttle protein (NSP) targets RLK family of genes to suppress antiviral response (Fontes et al. 2004). Thus, alteration in STK domains can confer resistance against viruses. The same is supported by another study, in which bamboo mosaic virus (BaMV) infection, upregulated the expression of STK gene of Nicotiana benthamiana and mutational knockdown of the same gene reduced BaMV titer in inoculated tissues (Cheng et al. 2013). Furthermore, differential transcriptome expression study conducted on resistant and susceptible blackgram genotypes showed upregulation of LRR receptor kinase, which led to strong conclusion regarding role of RLKs against viruses (Kundu et al. 2019; Raizada and Jegadeesan 2020). STKs presented plausible evidence of their role in viral defense, and the same was reflected by predominance of STK domains in qMYMIV6.1.1 suggesting the region as epitome of viral defense in blackgram genome.
General transcription factors such as Transcription initiation factor IIB (TFIIB) interacts with viral factors and regulates initiation of transcription (O’Brien and Ansari 2022). Therefore TFIIB is common target during pathogenesis of viruses and it was anticipated as functional target of virus in KUG253 and the same gene could be loss of function mutation in Mash114 because of presence of non-synonymous mutations in resistant parent. Alteration in sequence of TFIIB rendering it non-functional makes a strong potential candidate gene against viruses.
Some of the virus infections can induce host self-defense against invasion by stress response of endoplasmic reticulum (ER). During this process, there is active regulation of viral and host gene expression and switching of initiation of translation from cap dependent pathway to internal ribosome entry sites (IRES) dependent is often encountered. A cap dependent pathway initiates translation under normal conditions and IRES-dependent mechanism operates in certain conditions which include heat shock, infection and starvation, thus factors involved in IRES pathway could be strong candidates for resistance mechanism in Mash114.
The KASP markers developed from aforementioned 9 candidate genes were used for validation of QTL-seq results. KASP genotyping of RILs and parents identified a candidate region of 0.5 Mbp (1.9cM) with LOD 42.5 involved in virus resistance. The QTL region (1.9 cM) in genetic map spanned VM602, VM605 and VM610 markers originating from receptor protein serine/threonine kinase signaling pathway and defense responsive genes. VM602 and VM610 were closely placed in physical map at distance of only 28 kb, however the VM605 marker was present at distance of 0.55 Mbp from them. Furthermore, contrary to anticipation, some markers showed different order of placement on physical map in contrast to genetic map, thus implicating evolutionary cold-spots and hotspots for cross-overs across the chromosome length. The three closest markers i.e. VM602, VM605 and VM610 linked to MYMIV resistance gene carried 10 non synonymous SNP mutations between KUG253 and Mash114 in exonic regions of the potent candidate genes.
In nutshell, QTL-seq approach is very powerful technique for mapping desirable traits particularly in crops with limited genomic resources. The study has illuminated path for mapping traits of interest in blackgram for its genetic improvement and climate resilience. Serine/threonine kinases and plant NLR complexes are strong candidates against viral infection and they can be targeted as candidate host shields. Future strategy will be to dissect the putative candidate genes by plausible evidence through more recombination events in introgressed segment, functional validation and over-expression of the same by gene editing or transgenic approach and moving towards map based cloning of MYMIV resistance. Hence, our study has identified a major effect locus qMYMIV6.1.1 introgressed from ricebean and markers tightly linked to MYMIV resistance. This locus can be used by the breeders for developing blackgram varieties resistant to MYMIV using marker-assisted selection. Besides, given the crossing barriers between ricebean and other Vigna species, marker assisted transfer of introgressed segement into suitable genotypes of greengram and blackgram would be helpful for breeders using introgression line Mash114.
Funding
The funding of research was supported by Department of Biotechnology under project entitled as “Genetic Enhancement of Minor Pulses: Characterization, Evaluation, Genetic Enhancement and Generation of Genomic Resources for Accelerated Utilization and Improvement of Minor Pulses”.
Competing Interests
The author(s) declare that they have no competing interests.
Author contributions
SKD performed experiments and contributed to writing of manuscript, SK and RKG conceptualized the work and edited manuscript and finalized it, SK, AS and DB assisted in carrying out experiments, DB and AK contributed in manuscript editing.
Data Availability
All the original data is available in NCBI SRA database with accession number PRJNA788646
Acknowledgement
We acknowledge efforts of Dr. Inderjit Singh, Dr J. S. Sandhu, Dr. S. Gupta and Dr. Sarvjeet Singh for developing the genotype, ‘Mash114’.
- Afzal AJ, Wood AJ, Lightfoot DA (2008) Plant receptor-like serine threonine kinases: Roles in signaling and plant defense. Mol Plant-Microbe Interact 21:507–517 doi: 10.1094/MPMI-21-5-0507
- Arumuganathan K, Earle ED (1991) Nuclear DNA Content of Some Important Plant Species Nuclear DNA Content Material and Methods. Plant Mol Biol Report 9: 208–218
- Basak J, Kundagrami S, Ghose TK, Pal A (2005) Development of yellow mosaic virus (YMV) resistance linked DNA marker in Vigna mungo from populations segregating for YMV-reaction. Mol Breed 14:375–383 doi: 10.1007/s11032-005-0238-6
- Bashir M, Ahmad Z, Mansoor S (2006) Occurrence and distribution of viral diseases of mungbean and mashbean in Punjab, Pakistan. Pakistan J Bot 38:1341–1351
- Biswas KK, Malathi VG, Varma A (2008) Diagnosis of symptomless yellow mosaic begomovirus infection in pigeonpea by using cloned Mungbean yellow mosaic India virus as probe. J Plant Biochem Biotechnol 17:9–14 doi: 10.1007/bf03263253
- Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30:2114–2120 doi: 10.1093/bioinformatics/btu170
- Cao A et al. (2011) Serine/threonine kinase gene Stpk-V, a key member of powdery mildew resistance gene Pm21, confers powdery mildew resistance in wheat. Proc Natl Acad Sci U S A 108: 7727–7732 doi: 10.1073/pnas.1016981108
- Cheng SF et al. (2013) Ser/Thr Kinase-Like Protein of Nicotiana benthamiana Is Involved in the Cell-to-Cell Movement of Bamboo mosaic virus. PLoS One 8: 1–10 doi: 10.1371/journal.pone.0062907
- Churchill GA, Doerge RW (1994) Empirical threshold values for quantitative trait mapping. Genetics 138:963–971 doi: 10.1093/genetics/138.3.963
- Danecek P et al. (2011) The variant call format and VCFtools. Bioinformatics 27:2156–2158
- Fontes EPB, Santos AA, Luz DF, Waclawovsky AJ, Chory J (2004) The geminivirus nuclear shuttle protein is a virulence factor that suppresses transmembrane receptor kinase activity. Genes Dev 18:2545–2556 doi: 10.1093/genetics/138.3.963
- Gupta S, Gupta DS, Anjum TK, Pratap A, Kumar J (2013) Inheritance and molecular tagging of MYMIV resistance gene in blackgram (Vigna mungo L Hepper). Euphytica 193:27–37 doi: 10.1007/s10681-013-0884-4
- Hanson P J et al. (2012) IRES-dependent translational control during virus-induced endoplasmic reticulum stress and apoptosis. Front Microbiol 3:1–16 doi: 10.3389/fmicb.2012.00092
- Hurni S et al (2015) The maize disease resistance gene Htn1 against northern corn leaf blight encodes a wall-associated receptor-like kinase. Proc Natl Acad Sci U S A 112:8780–8785 doi: 10.1073/pnas.1502522112
- Illa-Berenguer E, Van Houten J, Huang Z , van der Knaap E (2015) Rapid and reliable identification of tomato fruit weight and locule number loci by QTL-seq. Theor Appl Genet 128:1329–1342 doi: 10.1007/s00122-015-2509-x
- Karthikeyan AS et al. (2004) Analysis of an isolate of Mungbean yellow mosaic virus (MYMV) with a highly variable DNA B component. Arch Virol 149:1643–1652 doi: 10.1007/s00705-004-0313-z
- Kumar S , Pandey G (2020) Biofortification of pulses and legumes to enhance nutrition. Heliyon 6:4–9 doi: 10.1016/j.heliyon.2020.e03682
- Kundu A, Singh PK, Dey A, Ganguli S, Pal A (2019) Complex molecular mechanisms underlying MYMIV-resistance in Vigna mungo revealed by comparative transcriptome profiling. Sci Rep 9 doi: 10.1038/s41598-019-45383-w
- Lee KJ, Kim K (2015) The rice serine/threonine protein kinase OsPBL1 (ORYZA SATIVA ARABIDOPSIS PBS1-LIKE 1) is potentially involved in resistance to rice stripe disease. Plant Growth Regul 77:67–75 doi: 10.1007/s10725-015-0036-z
- Lei L, Zheng H, Bi Y, Yang L, Liu H, Wang J, Sun J, Zhao H, Li X, Li J, Lai Y, Zou D (2020) Identification of a Major QTL and Candidate Gene Analysis of Salt Tolerance at the Bud Burst Stage in Rice (Oryza sativa L) Using QTL-Seq and RNA-Seq. Rice 13:1–14 doi: 10.1186/s12284-020-00416-1
- Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25:1754–1760 doi: 10.1093/bioinformatics/btp324
- Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R (2009) The Sequence Alignment/Map format and SAMtools. Bioinformatics 25:2078–2079 doi: 10.1093/bioinformatics/btp352
- Magwene PM, Willis JH, Kelly JK (2011) The statistics of bulk segregant analysis using next generation sequencing. PLoS Comput Biol 7:1–9 doi: 10.1371/journal.pcbi.1002255 doi: 10.1371/journal.pcbi.1002255
- Maiti S, Paul S, Pal A (2012) Isolation, characterization, and structure analysis of a non-TIR-NBS-LRR encoding candidate gene from MYMIV-resistant Vigna mungo. Mol Biotechnol 52:217–233 doi: 10.1007/s12033-011-9488-1
- Malathi VG and John P (2008) Gemiviruses infecting legumes. In: Govind P, Rao P, Kumar PL, Holguin-Pena RJ (ed) Characterization, Diagnosis & Management of Plant Viruses 97-123. Stadium Press LLC, Texas, USA.
- Mansfeld BN, Grumet R (2018) QTLseqr: An R Package for Bulk Segregant Analysis with Next‐Generation Sequencing. Plant Genome 11:180006 doi: 10.3835/plantgenome2018.01.0006
- McKenna A et al. (2010) The genome analysis toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20:1297–1303 doi: 10.1101/gr.107524.110
- Mishra G P et al. (2020) Yellow Mosaic Disease (YMD) of Mungbean (Vigna radiata (L) Wilczek): Current Status and Management Opportunities. Front Plant Sci 11:1–24 doi: 10.3389/fpls.2020.00918
- Muñoz-Amatriaín M et al. (2017) Genome resources for climate-resilient cowpea, an essential crop for food security. Plant J 89:1042–1054 doi: 10.1111/tpj.13404
- Naimuddin K, Akram M (2010) Detection of mixed infection of begomoviruses in cowpea and their molecular characterization based on CP gene sequences. Journal of Food Legume 23:191-95.
- Nariani, T. K. (1960) Yellow mosaic of mung (Phaseolus aureus). Indian Phytopathol 13: 24-29
- Singh, G. and Bhan, L. K. (1998) Diseases of mungbean and urdbean and their management. In: Upadhyay R (ed)., Pp. 311-71 IPM system in Agriculture
- O’Brien MJ, Ansari A (2022) Beyond the canonical role of TFIIB in eukaryotic transcription. Curr Genet 68:61–67 doi: 10.1007/s00294-021-01223-x
- Pandiyan M et al. (2008) Broadening the Genetic Base and Introgression of MYMV Resistance and Yield Improvement through Unexplored Genes from Wild Relatives in Mungbean. Plant Mutat Reports 2(1)
- Pootakham W et al. (2021) A chromosome-scale assembly of the black gram (Vigna mungo) genome Mol Ecol Resour 21:238–250 doi: 10.1111/1755-0998.13243
- Project Coordinators Report (2018) All India Coordinated Research Project on MULLaRP (Mungbean, Urdbean, Lentil, Lathyrus, Rajmash, Fieldpea) Kalyanpur, Kanpur: ICAR-Indian Institute of Pulses Research 1–46.
- Qazi J, Ilyas M, Mansoor S, Briddon RW (2007) Legume yellow mosaic viruses: Genetically isolated begomoviruses. Mol Plant Pathol 8:343–348 doi: 10.1111/j.1364-3703.2007.00402.x
- Raizada A, Jegadeesan S (2020) Comparative Transcriptomic Analysis Revealed Complex 2 Molecular Mechanisms Underlying Pests, Pathogens Resistance 3 and Seed Development in Wild and Cultivated Blackgram. bioRxiv 1109374041 doi: https://doiorg/101101/20201109374041
- Ramos A, Fu Y, Michael V, Meru G (2020) QTL-seq for identification of loci associated with resistance to Phytophthora crown rot in squash. Sci Rep 10:1–8 doi: 10.1111/pbi.13266
- Rogers SO, Bendich AJ (1989) Extraction of DNA from plant tissues. Plant Mol Biol Man 10:73–83 doi: 10.1007/978-94-009-0951-9_6
- Shen X et al. (2020) The NLR Protein Encoded by the Resistance Gene Ty-2 Is Triggered by the Replication-Associated Protein Rep/C1 of Tomato Yellow Leaf Curl Virus. Front Plant Sci 11:545306 doi: 10.3389/fpls.2020.545306
- Singh I, Sandhu JS, Gupta SK, Singh S (2013) Introgression of productivity and other desirable traits from ricebean (Vigna umbellata) into black gram (Vigna mungo). Plant Breed 132:401–406 doi: 10.1111/pbr.12068
- Souframanien J , Gopalakrishna T (2006) ISSR and SCAR markers linked to the mungbean yellow mosaic virus (MYMV) resistance gene in blackgram [Vigna mungo (L) Hepper]. Plant Breed 125:619–622 doi: 10.1111/j.1439-0523.2006.01260.x
- Stanke M et al. (2006) AUGUSTUS: A b initio prediction of alternative transcripts. Nucleic Acids Res 34:435–439 doi: 10.1093/nar/gkl200
- Takagi H et al. (2013) QTL-seq: Rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations. Plant J 74:174–183 doi: 10.1111/tpj.12105
- Tuberosa R , Salvi S (2007) From QTLS to Genes Controlling Root Traits in Maize. Scale Complex Plant Syst Res 15–24 doi:101007/1-4020-5906-x_2
- Varma A , Malathi VG (2003) Emerging geminivirus problems: A serious threat to crop production. Ann Appl Biol 142:145–164 doi: 10.1111/j.1744-7348.2003.tb00240.x
- Wang D et al. (2021) Bulked QTL-Seq identified a major QTL for the awnless trait in spring wheat cultivars in Qinghai, China. Biotechnol Biotechnol Equip 35:124–130 doi: 10.1080/13102818.2020.1857661
Table 1
Whole genome re-sequencing data statistics of parents and bulks
|
Sample |
Raw reads (M reads) |
Raw data (Gb) |
Clean Reads (M reads) |
Clean Data (Gb) |
Average depth (x) |
|
Mash114 |
40 |
6.0 |
32.6 |
4.9 |
10.5 |
|
KUG253 |
29.3 |
4.4 |
24.3 |
3.6 |
7.7 |
|
Resistant bulk (MYMIV-R) |
87.3 |
13.1 |
73.8 |
10.9 |
22.9 |
|
Susceptible bulk (MYMIV-S) |
96 |
14.4 |
82.2 |
11.4 |
25.2 |
Table 2
Quantitative trait loci (QTL) regions associated with resistance to mungbean yellow mosaic India virus in Mash114 identified through QTL-seq
