Mapping population consisting of 162 RILs derived from an intraspecific cross of HPKM249 × HPK4 was evaluated for morphological and phonological traits. The mapping population was developed through Single Seed Descent (SSD) method from F₂ to F₈ generations. The parental lines (HPKM249 & HPK4) show contrasting characteristics to each other for the traits under study (Figure 1 and Table 1).

During the crop season, F₈ progenies were planted in the experimental areas of department of Agriculture Biotechnology, CSK HPKV, Palampur, H.P. (latitude, 32.11 and longitude, 76.53) and at Hill Agricultural Research & Extension Centre, Bajaura, H.P. (latitude, 31.85 and longitude, 77.16) using Augmented Block Design (ABD) with two replications and four checks, VLG-1, HPKM249, HPK4 and HPK317 repeated after 20 rows. Each 1-m line consisted of 10 plants spaced 30 cms apart. Standard agronomic practices were followed for raising healthy crop. Phenotyping of various phenological and morphological traits were carried out at these locations (Palampur and Bajaura) over a period of two seasons (2016 and 2017). Data on five randomly taken plants in between the first and last plants of the same line was recorded. Three phenological traits viz. days to 50% flowering (FL), days to maturity (MT) and reproductive period (RP) and three morphological traits viz. plant height (PH, cm), number of primary branches (PB), and number of secondary branches (SB) were recorded. Plant height, number of primary branches, and number of secondary branches were measured just before the physiological maturity of the plant by taking readings on five plants in each line and averaging before analysis. Days from the date of sowing to the date when 50% of the plants in a line showed the first fully open flower and days from sowing to physiological maturity when 90% of the plants had turned brown were recorded for calculating days to 50% flowering (FL) and days to maturity (MT) respectively. The reproductive growth period (RP) was calculated as the days between the start of flowering and physiological maturity.

Traits distribution was studied using skewness and kurtosis statistics by Past 3.25 software. Pearson correlation coefficients and frequency distribution among different traits were calculated using the same software. The parental phenotypic variance was analyzed using ANOVA.

Young leaf tissues (0.5–1 g) of 162 RILs individuals along with parents (HPK4 and HPKM249) were used for isolation of genomic DNA using the modified cTAB method (Murray and Thompson, 1980). Concentration (ng/µl) and purity of isolated DNA was checked on agarose gel electrophoresis and quantified on a microvolume spectrophotometer (Biospec-nano, Shimadzu Biotech, United States) using Tris EDTA as blank. PCR-based markers from different sources were used to screen polymorphism between parent HPK4 and HPKM249. A total of 2011 PCR primers consisting of 63 Expressed Sequence Tag Simple Sequence Repeats (EST SSRs), 403 genic Simple Sequence Repeats (gSSRs), 387 genomic Simple Sequence Repeats (geSSRs), 24 drought specific Simple Sequence Repeats (dsSSRs), 300 Simple Sequence Repeats (SSR) from other legumes, 450 Random Amplified Polymorphic DNA (RAPD) markers, and 384 Conserved Ortholog Set (COS) markers of Cook’s Lab UCD, United States were employed in the present study. The polymorphic primer pairs were then used for genotyping of the mapping population.

For amplification using SSR markers, a total reaction mixture of 10.0 µl volume was prepared to contain 4.80 µl of sterilized distilled water, 2.0 µl of template DNA (13 ng/μl), 0.5 µl each of forward and reverse primer (5 µM), 0.5 µl of MgCl2 (25 mM), 1.0 µl of 10X PCR buffer (10 mM Tris-HCl, 50 mM KCl, pH 8.3), 0.5 µl of dNTP mix (0.2 mM each of dATP, dGTP, dCTP and dTTP) and 0.2 µl of Taq polymerase (5U/µl). The reaction mixture for RAPD markers consisted of 1.5 of µl DNA, 2.5 µl of 10X PCR buffer (10 mM Tris, pH 9.0, 50 mM KCl, 0.01% Gelatin, 1.5 mM MgCl2), 0.2 µl of dNTPs (25 mM), 1.0 µl of MgCl2 (25 mM), 1.0 µl of primer (2 μM/μl), 0.1 µl of Taq DNA polymerase (5 U/µl) and made final volume of 25 µl by adding deionised water. The PCR amplification of genomic DNA was performed on Veriti 384^® Thermal Cycler (Applied Biosystems, CA, United States). The PCR program for amplification was an initial denaturation cycle at 94°C for 5 min, followed by 45 cycles of denaturation at 94°C for 1 min, annealing at 50°C–65°C (for SSRs), and 37°C (for RAPDs) for 1 min and extension at 72°C for 2 min and final extension at 72°C for 7 min. The amplified products were resolved on either 6% Polyacrylamide Gel Electrophoresis in 1X TBE or 3% metaphor agarose gel (Lonza) in 1X TAE buffer depending on the size difference between amplified DNA along with the size markers (100 bp DNA ladder). The fragments were visualized using Gel-Documentation Unit (ENDUROTM GDS Gel Documentation System, United States) or silver-staining procedure.

The input file for linkage map construction was prepared manually by the scoring of amplified bands. HPKM249 type banding pattern was scored as “A”, HPK4 type banding pattern was scored as “B” and heterozygous loci was scored as “H”, whereas for RAPD markers the absence (0) and presence 1) of bands were recorded. Sizes of amplified fragments were noted by using a 100-bp DNA ladder (Fermentas, Lithuania).

The genetic linkage map construction was performed with scored genotypic data file using JOINMAP^® 4.1 program (van Ooijen, 2006). The LOD threshold >3.0 and <8.0 with a step-up of 0.5 was considered significant for identifying different linkage groups and clustering of markers on them. The linkage groups which show the highest number of markers with maximum linkage among them at different LODs values were selected.

For quantitative trait analysis, QTLs were identified using the composite interval mapping (CIM) method (Zeng, 1993; Zeng, 1994) implemented in Windows QTL Cartographer V2.5 software (Wang et al., 2005). The walking speed selected for QTLs was 2 cM with a window size of 10 cM using the Zmapqtl standard model 6. The forward regression algorithm was used to obtain cofactors. To calculate a genome-wide threshold for LOD score a 1000-permutation test at a significance level of p = 0.05 for shuffling genotypes with the phenotype means was performed (Doerge and Churchill, 1996). A LOD threshold score of ≥2.5 was selected for identification of the QTLs on the horsegram LGs. The location of QTLs with 95% confidence intervals was estimated by one LOD interval around the QTL peak (Mangin et al., 1994). The estimated additive effect and the percentage of phenotypic variation explained (based on the R2 value) by each putative QTL was calculated by the Zmapqtl procedure using the software with the CIM model. The QTL map was prepared with MapChart 2.32 software (Voorrips, 2002).

Horsegram being an orphan crop, very little has been understood about the genetic structure of its traits. We analyzed six phenological and morphological traits as described in Table 2 which represents the descriptive statistics for all the traits. Parental lines, HPKM249 and HPK4 which were used for developing the mapping population found to have significant differences for all the traits studied. Three phenological traits, namely days to 50% flowering (FL), reproductive period (RP), and days to maturity (MT) are important indicators of maturity and were used for phenotyping of the RILs population. Phenotyping of FL showed significant genetic variability for RILs in different years and locations. HPKM249 flowered in 36 days as compared to 54 days of HPK4 during 2016 at Palampur and similar results were observed at Palampur during 2017, whereas at Bajaura during 2017, HPKM249 flowered in 32 days as compared to 57 days of HPK4. The range for days to flowering among RILs varied from 30–58 days in 2016 at Palampur, 32–52 days in 2017 at Palampur, and 31–57 days in 2017 at the Bajaura location. Further, no significant difference was found among RILs in different years and locations. A similar trend was observed for RP and MT with no significant difference among RILs in different years and locations.

RIL population was phenotyped for various morphological traits like plant height (PH), number of primary branches (PB), and number of secondary branches (SB). Plant height (PH) varied from 34–98 cm and 48–106 cm at Palmapur for the year 2016 & 2017, respectively. A significant difference was observed for PH in 2017 at Bajaura (60–145 cm) as compared to the Palampur location. Similarly, PB varied from one to six branches at Palampur in both years (2016 and 2017), however a significant difference was observed for PB in 2017 at Bajaura, which varied from 3–10. Further, significant differences for PH and PB among RILs were observed in both seasons and locations. A similar result was observed for SB with significant differences among RILs in both years and locations.

The phenotypic values showed a normal frequency distribution, which is a typical characteristic of quantitative traits (Figure 2). Correlations values among traits in RIL populations showed that the traits were positively correlated with each other and were statistically significant (p < 0.05). The maximum correlation was found between MT and RP (r = 0.86) followed by PB and SB (r = 0.85) (Figure 3 and Table 3). All traits showed positive correlations among themselves except for FL and RP which showed negative correlation with each other in all environments. The ANOVA of 162 RILs for different locations and environments showed significant variation for all the traits.

The polymorphism survey was done on parental lines using 1177 SSR primer pairs [63 (Horsegram EST SSRs) + 403 (Horsegram genic SSRs) + 387 (Horsegram genomic SSRs) + 24 (drought specific SSRs) + 300 (SSRs from other legumes viz. red clover and Medicago)] of which 430 were found to be polymorphic (Table 4). Along with this, 450 RAPD primers were also screened out of which 55 were found to be polymorphic and in addition 384 COS markers of Medicago truncatula were screened out of which 8 were found to be polymorphic (Table 4). A total of 2011 primers were screened, of which 493 were found to be polymorphic and were further used for genotyping of 162 individuals (Table 4). The size of amplified bands produced by 493 polymorphic primers ranged between 100 and 250 bp. The genotyping data thus obtained was scored manually and used as the input file for the construction of the horsegram linkage map.

Using JoinMap software, version 4.0, of the total 493 polymorphic markers, 295 (59.84%) were mapped on 10 LGs at LOD 3.5. The linkage map spanned 1,541.7 cM (Kosambi cM) length and an average marker interval size of 5.20 cM (Table 5). These 295 mapped markers include 15 EST SSRs, 87 genomic SSRs, 87 genic SSRs, 22 RAPDs, 73 SSRs from other species, four drought-specific markers, and seven COS markers. Each of the ten linkage groups differed from one another in terms of their length and the total number of markers mapped. Of the total 295 mapped markers, LG1 harbored 89 markers followed by LG2 which contained 58 markers. 35 markers were mapped on LG3, 29 markers were mapped on LG4, 19 were mapped on LG5 and LG7 whereas 18 markers were mapped on LG6. LG eight and LG nine contained the least number of markers with seven and 6 markers, respectively, and linkage group 10 harbored 15 markers (Figure 4 and Table 5). Though LG7 is having a maximum length of 238.5 cM due to very less number of markers 19) present on it with an average marker density of 12.5 however is of less importance while mapping of different QTLs in comparison to linkage groups LG1 and LG2 having marker density of 2.0 and 2.7 respectively even having smaller linkage size (182.9 and 159 cM). Segregation distortion for all the 493 polymorphic markers was determined and of these 295 (59.83%) followed the expected segregation ratio, whereas 198 markers (40.16%) were found to show deviation. These distorted markers were excluded from the final analysis.

The maximum distance recorded between two markers was 61.2 cM on LG7 and the minimum distance was 0.003 cM on LG2. The number of markers present in different linkage groups was unequal. Four large groups having 12–19 markers within a length of 10 cM and five groups having 28–31 markers within a length of 30 cM were found. The length of the linkage groups did not reflect the number of markers linked on it as the distance between markers varied across different linkage groups. For example, LG1 carrying 89 markers having a length of 182.9 cM and an average marker distance of 2.0 cM, whereas LG4 having a length of 188.0 cM carrying 29 markers and an average marker distance of 6.5 cM, and LG5 having a length of 192.5 cM covered by 19 markers and average spacing of 10.0 cM.

The QTL analysis identified four QTLs (LOD ≥2.5) for phenological and seven QTLs (LOD ≥2.5) for morphological traits (Figure 4 and Table 6). One QTL for days to 50% flowering (qFT01) was detected on LG2 flanked by markers MUGR644-MUMST80 at LOD score of 2.8 and explaining 6.62% of the phenotypic variation with allelic contribution by HPKM249 resulted in reduced flowering time by about >2 days. One QTL for the reproductive period (qRP01) was detected on LG5 flanked by MUGSSR10-RCS6448 at a LOD score of 2.7 and explaining 6.36% of the phenotypic variation (Table 6). Two QTLs namely qMT01 and qMT02 were detected for days to maturity on LG7 flanked by markers MUGSSR241-HUGMS39 at LOD 2.6 explaining 7.25% of the phenotypic variation and on LG9 flanked by HUGMS3-MUGR607 marker interval at LOD 2.9 explaining 47.53% of the phenotypic variation, respectively and in combination explained 54.78% of total phenotypic variation. Additive effect is the difference in the average performance of the RIL carrying early maturity allele of first parent with respect to those carrying the allele of second parent at the particular locus (QTL). A positive value (+) of the additive effect indicates the allele originating from HPKM249 and a negative value (–) of the additive effect indicates the allele originating from HPK4.

Additive effect demonstrated that HPKM249 contributed alleles for a reproductive period and days to maturity with QTL named qMT01 resulted in reduced days to maturity by > 3days and qMT02 resulted in reduced days to maturity by > 8 days (Table 6).

One QTL was detected for plant height namely qPH01 on LG1 flanked by RCS6168-RCS6169 and explaining 6.6% of the phenotypic variation at a LOD value of 2.7. This QTL had an additive effect of 3.96 cm and was contributed by the allele from HPKM249. Three QTLs were detected for primary branches namely qPB01, qPB02, and qPB03 with two on LG6 (qPB01 and qPB02) both flanked by OPI66-MUMST29 and explaining 22.0 and 17.0% of the phenotypic variation, respectively, and one QTL on LG9 (qPB03) flanked by HUGMS3-MUGR607 explaining 32.4%of the phenotypic variation (Table 6). These QTLs in combination explained 71.4% of total phenotypic variation for primary branches. Additive effect demonstrated allelic contribution from both the parents. Three QTLs were detected for secondary branches all on LG7 namely qSB01, qSB02, and qSB03. Both qSB01 and qSB03 were flanked by MUD77-HUGMS13 marker interval explaining 23.6 and 15.5% phenotypic variation respectively. qSB02 was flanked by MUMS13-MUMS95 marker interval explaining 7.5% of phenotypic variation. These QTLs in combination explained 46.6% of total phenotypic variation for secondary branches. Additive effect demonstrated allelic contribution from both the parents. The position of QTLs for phenological and morphological traits at LOD ≥2.5 on 10 linkage groups on horsegram has also been shown in Figure 5.