High-Density GBS-Based Genetic Linkage Map Construction and QTL Identification Associated With Yellow Mosaic Disease Resistance in Bitter Gourd (Momordica charantia L.)

Yellow mosaic disease (YMD) in bitter gourd (Momordica charantia) is a devastating disease that seriously affects its yield. Although there is currently no effective method to control the disease, breeding of resistant varieties is the most effective and economic option. Moreover, quantitative trait locus (QTL) associated with resistance to YMD has not yet been reported. With the objective of mapping YMD resistance in bitter gourd, the susceptible parent “Punjab-14” and the resistant parent “PAUBG-6” were crossed to obtain F4 mapping population comprising 101 individuals. In the present study, the genotyping by sequencing (GBS) approach was used to develop the genetic linkage map. The map contained 3,144 single nucleotide polymorphism (SNP) markers, consisted of 15 linkage groups, and it spanned 2415.2 cM with an average marker distance of 0.7 cM. By adopting the artificial and field inoculation techniques, F4:5 individuals were phenotyped for disease resistance in Nethouse (2019), Rainy (2019), and Spring season (2020). The QTL analysis using the genetic map and phenotyping data identified three QTLs qYMD.pau_3.1, qYMD.pau_4.1, and qYMD.pau_5.1 on chromosome 3, 4, and 5 respectively. Among these, qYMD.pau_3.1, qYMD.pau_4.1 QTLs were identified during the rainy season, explaining the 13.5 and 21.6% phenotypic variance respectively, whereas, during the spring season, qYMD.pau_4.1 and qYMD.pau_5.1 QTLs were observed with 17.5 and 22.1% phenotypic variance respectively. Only one QTL qYMD.pau_5.1 was identified for disease resistance under nethouse conditions with 15.6% phenotypic variance. To our knowledge, this is the first report on the identification of QTLs associated with YMD resistance in bitter gourd using SNP markers. The information generated in this study is very useful in the future for fine-mapping and marker-assisted selection for disease resistance.


INTRODUCTION
Bitter gourd (Momordica charantia L.) is also popularly known as bitter melon, bitter cucumber, balsam pear, karela, and African pear (Basch et al., 2003). It belongs to the family Cucurbitaceae and genus Momordica having chromosome number 2x = 2n = 22 and is one of the important vegetable crops grown in the tropics. There is uncertainty in the exact place of the origin of bitter gourd but according to various molecular studies, east India is supposed as its center of origin, whereas China is its secondary center of diversity (Dey et al., 2006;Gaikwad et al., 2008). Among various bitter gourd growing regions, India, Malaysia, China, South America, and Africa are on the top of the list (Miniraj et al., 1993). This crop is mainly cultivated in warm climates, so the summer season is the most appropriate season for its production in India. Though all parts of the plant are edible, the immature fruits are consumed most commonly. Bitter gourd is rich in vitamins and nutrients (Raman and Lau, 1996;Blum et al., 2012) and possess antioxidant (Krishnaiah et al., 2011) and antihepatotoxic (Semiz and Sen, 2007) and anti-cancer properties (Brennan et al., 2012).
Due to these medicinal properties, people have shown interest in including M. charantia in their diet, and hence the demand for bitter gourd has increased (Warrier et al., 1997). But its production is decreasing due to various diseases caused by different pathogens like bacteria, fungi, and viruses. Among these, yellow mosaic disease (YMD) has proven to be the most destructive viral disease which could result in up to 100% yield loss depending upon the developmental stage of the crop (Giri and Mishra, 1986). Various efforts were made to study the transmission, causal pathogen, host range, and etiological impact of the disease. YMD is caused by bipartite geminivirus belonging to genus Begomovirus and family Geminiviridae, having circular topology and a unique single-stranded DNA genome. This group of viruses can easily manipulate their genome which helps them to infect a wide range of dicotyledonous plants. Begomoviruses also cause YMD in some other crops like cotton, mungbean, soybean, cucumber, okra, etc. The viruses affect the host at all growth stages and are transmitted by Bemisia tabaci (whitefly). The disease symptoms include typical mosaic, upward leaf curling, mottling, crinkling, and severe stunting symptoms (Matthew et al., 1991).
To date, many viruses related to YMD have been identified and characterized. These include Papaya ring spot virus (Gonzalez et al., 2003), Zucchini yellow mosaic virus (ZYMV) (Fukumotu et al., 1993), Water melon silver mottle virus (Tokashiki and Yasuda, 1991), Water melon mosaic virus-1 (Tomar and Jitendra, 2005), Bitter gourd yellow mosaic virus (Rajinimala et al., 2005), Tomato leaf curl distortion virus (Tahir and Haider, 2005), Cucumber mosaic virus (Takami et al., 2006), and Indian cassava mosaic virus (Rajinimala and Rabindran, 2007). Recently, Squash vein yellowing virus and Cucurbit leaf crumple virus (Adkin et al., 2008), Melon yellow spot virus (Takeuchi et al., 2009) and Pepper leaf curl Bangladesh virus  have also been investigated on bitter gourd. Association of Tomato leaf curl New Delhi virus with YMD of bitter gourd from India has been reported in isolation (Tiwari et al., 2010) or as a mixed infection with ZYMV by Sharma et al. (2015) while Ali et al. (2010) reported the presence of Tomato leaf curl Palampur virus in YMD affected bitter gourd plants. The simultaneous occurrence of different viruses in bitter gourd plants results in the mosaic complex. The ability of genome manipulation by recombination coupled with vector transmission of viruses by vectors to new areas has doubled the impact of the diseases (Padidam et al., 1999).
The chemical methods used to control the vector whitefly population generally pollute the environment and cause detrimental health issues. One of the most sustainable ways of managing the YMD in bitter gourd is through host plant resistance. So the development of improved disease-resistant cultivars should be the primary goal for the breeders but the major problem in YMD resistant breeding is the screening for the disease. The development of non-uniform disease due to the fluctuation in vector population is generally observed in different locations and seasons. This problem can be solved by Marker-assisted selection (MAS), as markers linked to traits of interest help to reduce the number, time, and cost of phenotypic evaluations. So far, there is no report on mapping and DNA markers associated with YMD resistance in bitter gourd.
An enriched genetic linkage map is essential to identify the location of genes controlling disease resistance. A few linkage maps have been constructed in bitter gourd by using DNA molecular markers. Kole et al. (2012) developed the first linkage map using amplified fragment length polymorphism (AFLP) markers. Simple sequence repeats (SSR), AFLP, and sequence-related amplified polymorphism (SRAP) markers were used in the linkage map developed by Wang and Xiang (2013). Owing to the development made in next-generation sequencing technologies, sequencing-based genotyping methods like Restriction enzyme-associated DNA sequencing (RAD-seq; Baird et al., 2008) and genotyping by sequencing (GBS; Elshire et al., 2011) can be used in single nucleotide polymorphism (SNP) markers discovery, linkage map development and quantitative trait locus (QTL) identification in bitter gourd. Urasaki et al. (2017) published the first short-read scaffold-level assembly of bitter gourd inbred line, OHB3-1, and constructed a linkage map consisting of 11 linkage groups using RAD-seq analysis. Cui et al. (2018) also used RAD-seq SNP markers for linkage map construction in bitter gourd and GBS without reference genome was used by Rao et al. (2018). Recently, a chromosome-level genome assembly is published along with a linkage map, and it is currently one of the most complete assemblies among publicly available genomes in cucurbitaceae (Matsumura et al., 2020).
In this study, we used a reference-based GBS approach for SNP identification and genotyping of a mapping population, segregating for YMD. The identified SNPs were used in the construction of linkage maps and detection of loci associated with YMD.

Plant Material
Resistant parent, PAUBG-6 and susceptible parent Punjab-14 were used for the development of the mapping population. The Department of Vegetable Science, Punjab Agricultural University (PAU), screened bitter gourd germplasm against YMD following field screening technique and identified bitter gourd genotypes with YMD resistance. Of these PAUBG-6 was found to be a stable source of resistance to this viral disease. Whereas the recipient Punjab-14 is an elite cultivar with green shining fruits suitable for cooking through stuffing. Female flowers of Punjab-14 were hand-pollinated with pollen from male flowers of PAUBG-6 to generate F 1 plants. F 1 s were further hand self-pollinated to produce F 2 seeds and individual F 2 plants were planted again. A total of 101 F 2 plants were selfed individually to produce F 2:3 lines. These lines were selfed in subsequent generations to develop F 2:4 and F 2:5 lines (Supplementary Figure 1). The populations were grown according to the standard agronomic practices in the spring and rainy seasons of the years 2017-2020.

Phenotypic Screening and Disease Evaluation
F 4:5 individuals were phenotyped for YMD resistance. The F 4:5 individuals were screened using both artificial and natural inoculations to achieve a uniform vector distribution. A total of 101 F 4:5 families were subdivided into three subsets, which were then grown in three different environments consisting of three replications in randomized block design. The first subset was grown in a nethouse (2019), the second during the rainy season (2019), and the third during the spring season (2020) to study the disease response. For artificial inoculation, the nonviruliferous whiteflies (B. tabaci) were reared on virus-free cotton plants in an insect-proof cage of fine mesh held in a controlled temperature chamber at 28 • C with a photoperiod of 14 h and 30-50% relative humidity (Polston and Capobianco, 2013). The virus-infected bitter gourd plants showing YMD symptoms were collected from the field and were used as a primary source of the virus. The non-viruliferous whiteflies were collected from cotton plants and fed on virus-infected bitter gourd plants for virus acquisition for 24 h. Later these viruliferous whiteflies were released on test plants sown in plug trays and maintained in the same cage. The F 4:5 inoculated individuals were then transplanted in testing environments for screening. The spreader row technique was adopted in screening to evaluate lines against YMD infection. About 10 F 5 seedlings from each F 4 plant were planted on well-prepared hills of 1.5 m wide beds with a plantto-plant distance of 45 cm. Punjab-14 was planted as indicatorcum-infestor row after every four test rows. The main source of virus transmission is the whitefly; to maintain the natural population of whiteflies in the experimental field, insecticide spray was avoided. In field conditions during the rainy and spring seasons, plants were under constant disease pressure while the plants in the nethouse were exposed to viruliferous whiteflies only once before transplanting. The disease reaction was recorded at the time of the first harvest using a grading scale suggested by Arunachalam et al. (2002). The lines with a score of 0 were considered as highly resistant, a score of 1 as resistant, a score of 2 as moderately resistant, a score of 3 as moderately susceptible, a score of 4 as susceptible, and a score of 5 as highly susceptible (Supplementary Figure 2). The percent disease index (PDI) was calculated for each line by applying the formula given by McKinney (1923): Here, s: disease grade; n: number of plants in the disease grade, S: highest disease rating scale, and N: total numbers of plants.

DNA Isolation, Library Construction, and Sequencing
Genomic DNA was isolated from fresh leaf tissue of F 4 plants as well as from parents for genotyping by using CTAB (Cetyl Trimethyl Ammonium Bromide) method as given by Doyle and Doyle (1990). The DNA was quantified using NanoDrop TM 1000 spectrophotometer (Thermo Scientific, Wilmington, United States) whereas, the integrity and quality of DNA was checked on 0.8% agarose gel. High-quality DNA diluted to ∼100 ng/µl was used for GBS analysis. DNA digestion, adapter ligation, library construction, and sequencing were carried out by Novogene (China). For bitter gourd library preparation, two restriction enzymes Mse1 and Msp1 were used. The library was sequenced with the paired-end (PE) read length of 150 bp using Illumina HiSeq 2000.

Single Nucleotide Polymorphism Identification and Genotyping
The PE sequence reads from FASTQ files were filtered based on criteria for sequence quality and length using a trimmomatic tool (Bolger et al., 2014). Sequence reads having < 50 bp of length were removed. The filtered reads were aligned to the reference genome of genotype sequence OHB3-1 of bitter gourd using Bowtie 2 software (Langmead and Salzberg, 2012) with default parameters. Sequence alignment map (SAM) to binary alignment map (BAM) format conversion followed by sorting was done using SAMtools. SNP calling was performed using the BCFtools "mpileup" command. The resulting SNPs were filtered out with 25% missing data, minimum read depth 4, minimum base quality 20, and minor allele frequency (MAF) of ≤0.05. Then, further filters like SNPs with any missing site in one or both parents and any site with a heterozygous call at one or both parents were also removed. The monomorphic SNP sites were rejected and only polymorphic SNP sites called in both the parents were selected. A total of 4196 SNPs with contrasting alleles were retained after converting the generated SNPs into a parent-based format and were further used for linkage map construction.

Genetic Linkage Map Construction
The filtered SNPs from the 101 lines were used for the construction of the linkage map using Mapdisto version 1.7.7 (Lorieux, 2012). The chi-square test (χ 2 ) was performed to check the segregation of each SNP marker for any deviation from the expected Mendelian ratio. Highly distorted and unlinked markers were excluded from the analysis. The SNPs were split into linkage groups with a recombination frequency of 0.3 at a logarithm of odds (LOD) value of 3-6. The markers were ordered and rippled by the Seriation method and SARF (Sum of Adjacent  Recombination Frequencies) criteria respectively. Finally, the Kosambi mapping function was applied to compute map distance (cM). Linkage groups (LGs) were drawn using MapChart V2.2 (Voorrips, 2002).

Quantitative Trait Locus Analysis
The QTL analysis was conducted using the IciMapping software program (Meng et al., 2015) by ICIM-ADD (Inclusive composite interval mapping) method. Only QTLs having LOD values > 3.0 and phenotypic variance (PV) of more than 10% were considered putative QTLs for YMD resistance. Nomenclature for the QTL was denoted as, for instance, qYMD_pau4.1 where "qYMD" represents the QTL for resistance to YMD, "pau" represents institution name, and 4_1 represents the first QTL on chromosome 4.

Phenotypic Evaluation of Mapping Population
The disease development was very good in three environments. The YMD symptoms like stunted growth, shortening of internodes, deformed fruits, and increase in hairiness were observed in the experimental field (Figure 1). The parents of the mapping population showed contrasting phenotypic reactions for YMD resistance. In F 1 all plants were susceptible indicating disease resistance to be recessive. In the case of the F 4:5 population, the mean PDI score ranged from 6 to 82, 9 to 98, and 6 to 98 with a coefficient of variation (CV) of 42.1, 38.9, and 44.7% for the nethouse, rainy, and spring seasons respectively ( Table 1). The frequency distributions (Figure 2) and the absolute value of <1 for skewness and kurtosis (Table 1) for the trait demonstrate the continuous and normal distribution of the disease reaction data.

High-Throughput Sequencing and SNP Identification
Sequencing of 101 F 4 lines generated over 144 million PE reads (21.73 GB) that had a length of 150 bp. In parents, a total of 1.8 and 1.5 million reads in Punjab-14 and PAUBG-6 were generated respectively. The reads from F 4 individuals ranged from 0.5 to 3.3 million reads with an average of 1.4 million reads (Supplementary Figure 3). Then, the cleaned 88-99% sequencing reads were successfully aligned with a reference genome (DNA Data Bank of Japan, accession nos. BLBB01000001-BLBB01000011) for SNP identification. As a result, 4,31,148 raw SNPs were identified, across the 11 chromosomes. Out of these, 1, 54,552 SNPs were polymorphic with a polymorphism rate of 35.84%. The 17,488 SNPs were filtered with filtering 5% MAF and 25% missing data. These SNPs were then filtered for Chi-Square at p = 0.001 which resulted in 4196 SNPs ( Table 2). These 4196 SNPs were then used for linkage map construction. The distribution of SNPs on 11 chromosomes is illustrated in Figure 3.

Single Nucleotide Polymorphism Based High-Density Linkage Maps
Out of 4,196 SNPs, 3,144 formed 15 linkage groups spanning a cumulative length of 2415.2 cM (Figure 4). Linkage groups were named according to chromosomes. The number of linkage groups obtained was more than the haploid number of chromosomes in the bitter gourd genome. This was due   Table 3).

Detection of QTLs Conferring YMD Resistance
The linkage map developed from genotyping of 3,144 SNP markers was integrated with phenotypic disease reaction data from three independent trails and was analyzed using the ICIM method. Three QTLs for YMD resistance were identified on chromosome 3, 4, and 5 and designated as qYMD_pau3.1, qYMD_pau.4.1, and qYMD_pau.5.1 respectively ( Table 4). The QTL qYMD_pau3.1 was identified for the rainy season at position 149 cM between flanking SNP markers C3_11081069 and C3_14788898, showing an LOD score of 3.05 (Figure 5). This QTL explained 13.5% phenotypic variance (PV) with an additive effect of 5.60. The QTL qYMD_pau.4.1 flanked by C4_8810605-C4_8985192 markers was detected in the rainy season at 104 cM with an LOD score of 4.97 and in the spring season with an LOD score of 4.29 on chromosome 4. This QTL showed an additive effect of −6.82 and accounted for 17.5-21.6% disease score variation for both seasons. The QTL qYMD_pau.5.1 was identified in the spring season with a high LOD score of 5.23 and 22.1% phenotypic variance at position 152 cM on chromosome 5 with flanking markers C5_7937924-C5_8770358 whereas for the nethouse QTL was detected on chromosome 5 at position 148 cM flanked by C5_8409196-C5_8409575 with LOD value of 3.87 and explained 15.6% phenotypic variance. This QTL lies within the same confidence interval as the QTL region identified in the spring season on chromosome 5.

DISCUSSION
Yellow mosaic is one of the severe threats in bitter gourd production that needs to be controlled by using a sustainable, environment-friendly approach. Breeding of resistant varieties through MAS is the most effective way as it has proved its importance in achieving desirable phenotypic variation in less time. There is no report of inheritance and mapping of YMD resistance gene in bitter gourd to date and hence no marker linked to the trait has been reported. Therefore, this study aimed at identifying YMD resistance genes through the GBS approach. GBS generates a large number of SNPs, which are used in the construction of high-density linkage maps and gene tagging. This approach has been used in cucurbits like watermelon (Branham et al., 2016;Meru and McGregor, 2016), bitter gourd (Rao et al., 2018), bottle gourd (Wu et al., 2017), squash , musk melon (Chang et al., 2017;Daley et al., 2017;Branham et al., 2018), and pumpkin (Montero-Pau et al., 2017) for construction of linkage map and QTL mapping. Biparental mapping in integration with SNP genotyping is becoming a powerful tool for dissecting complex traits in plants (Xu et al., 2017). In this regard, a genetic map with reasonable marker density and genetic regions associated with YMD resistance identified in the present study is an important milestone in bitter gourd breeding. We used GBS to discover and genotype SNPs in F 4 mapping populations, derived from the cross between Punjab-14 (susceptible parent) and PAUBG-6 (resistant parent). Sequencing of mapping population resulted in the identification of 144 million reads in the F 4 population with an average of 1.4 million reads per individual. The 88-99% alignment rate with the reference genome indicates the sequencing accuracy. We used chromosome level genome assembly (Matsumura et al., 2020) as a reference to find SNPs on bitter gourd chromosomes while in previous reports either the short reads scaffold assembly (Cui et al., 2018) or non-reference genome (Rao et al., 2018) approach was used for SNP discovery. About 4,31,148 raw SNPs were identified which were passed through stringent filter criteria like missing percentage, MAF, and segregation behavior to select SNPs for constructing high-density genetic maps. This resulted in the reduction of SNPs from millions to few thousands, however, this is very common in GBS data, and similar results have been reported in many other GBS-based genetic mapping studies Mathivathana et al., 2019). In summary, in the present study, a total of 3,144 SNPs were finally mapped on the linkage map. The number of mapped markers was low in comparison to the identified SNPs between parental lines  Frontiers in Plant Science | www.frontiersin.org may be due to the limited sequencing depth used for GBS. The linkage map spanning a total 2415.2 cM length consisted of 15 linkage groups with a mean marker interval of 0.7 cM an average of 209.6 markers for each linkage group. The mean marker density (marker/cM) 1.30 is higher than 0.86 reported by Rao et al. (2018), 0.30 by Matsumura et al. (2014), 0.42 by Urasaki et al. (2017), and 0.46 by Cui et al. (2018). The linkage map produced in this study resulted in more linkage groups than haploid chromosome number (n = 11). This may be due to the absence of linkage between SNPs in that region due to a relatively small F 4 population size. The population development in bitter gourd crops is laborious and time-consuming, as the crop is highly cross-pollinated and monoecious in nature. But efforts were made to advance the population from F 2 to F 2:4 and F 2:5 lines to achieve more homozygosity and accurate results. Infection at the early developmental stage results in 100% loss; to avoid this situation, phenotyping for the disease resistance was performed in the F 4:5 population. The virus incidence and disease severity vary with varying seasons and environments in India. The temperature during the spring season is generally high with low humidity as compared to the high humidity of the rainy season, which helps the vector population to grow at a faster rate. This results in increasing the whitefly population and virus incidence. Mapping populations in the present study were phenotyped in three independent environment conditions. In the nethouse, the plants were challenged with viruliferous whiteflies for once before transplanting, while during the rainy season (2019) and spring season (2020) the plants were under constant disease pressure. Generally, seasonal and annual variations in the field influence the virus incidence and severity of the disease, so to minimize these variations, artificial inoculations are done before transplanting the test material in fields for further screening (Picó et al., 1998). So to avoid the nonuniform development of disease while screening, artificial and natural whitefly mediated transmission techniques were used. The phenotypic data depicted the polygenic inheritance of the disease. Huge variations were observed among progenies across the test environments which can be seen in mean PDI scores ( Table 1). The phenotypic variation ultimately affects the identification of stable QTLs. This might be due to the presence of environment interactions/pathogenic variability across the seasons.
Quantitative trait locus mapping for YMD resistance has revealed three significant QTLs qYMD_pau3.1, qYMD_pau.4.1, and qYMD_pau.5.1 having major effects with PV of more than 10% on chromosomes 3, 4, and 5. The identified QTLs in the present study are novel as this is the first report of mapping QTLs for YMD resistance in bitter gourd. Out of these, QTL qYMD_pau.4.1 was identified on chromosome 4 in two seasons. The QTL qYMD_pau3.1 is observed during the rainy season, while qYMD_pau5.1 is observed during the spring season and under nethouse conditions. The additive effect of qYMD_pau4.1 and qYMD_pau5.1 indicates that the alleles for resistance are contributed by the resistant parent (PAUBG-6) while the alleles for the resistance due to qYMD_pau3.1 are contributed by the susceptible parent . Altogether, these results suggested that the dense genetic map provides accurate detection of QTLs for YMD resistance in bitter gourd. The inclusive composite interval mapping (ICIM-ADD) method was used in the identification of QTLs, as it has increased detection power, produces a lower false detection rate, and avoids the possible sampling variance and the complicated background marker selection process while retaining all the advantages of CIM.

CONCLUSION
In this study, SNPs were identified from the bitter gourd genome by using the GBS approach. Using the identified SNPs a linkage map was developed between Punjab-14 and PAUBG-6. Three QTLs for YMD resistance have been mapped successfully, and this is the first report of mapping YMD resistance in bitter gourd to the best of our knowledge. The information generated on QTLs is very useful for its fine-mapping and for MAS for the trait in the future. These can be applied for the development of YMD resistant genotypes and molecular dissection of YMD resistance in bitter gourd.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA702876.

AUTHOR CONTRIBUTIONS
GK performed the YMD phenotyping, contributed to the analysis of SNP data, linkage map construction, and QTL analysis. AS contributed to the YMD phenotyping and data analysis. DS contributed to the SNP identification. NS and MP designed and supervised the project. NS edited the manuscript. PC gave critical suggestions throughout the study and also edited the manuscript. All authors contributed to the article and approved the submitted version.