Molecular Mapping of Flowering Time Major Genes and QTLs in Chickpea (Cicer arietinum L.)

Flowering time is an important trait for adaptation and productivity of chickpea in the arid and the semi-arid environments. This study was conducted for molecular mapping of genes/quantitative trait loci (QTLs) controlling flowering time in chickpea using F2 populations derived from four crosses (ICCV 96029 × CDC Frontier, ICC 5810 × CDC Frontier, BGD 132 × CDC Frontier and ICC 16641 × CDC Frontier). Genetic studies revealed monogenic control of flowering time in the crosses ICCV 96029 × CDC Frontier, BGD 132 × CDC Frontier and ICC 16641 × CDC Frontier, while digenic control with complementary gene action in ICC 5810 × CDC Frontier. The intraspecific genetic maps developed from these crosses consisted 75, 75, 68 and 67 markers spanning 248.8 cM, 331.4 cM, 311.1 cM and 385.1 cM, respectively. A consensus map spanning 363.8 cM with 109 loci was constructed by integrating four genetic maps. Major QTLs corresponding to flowering time genes efl-1 from ICCV 96029, efl-3 from BGD 132 and efl-4 from ICC 16641 were mapped on CaLG04, CaLG08 and CaLG06, respectively. The QTLs and linked markers identified in this study can be used in marker-assisted breeding for developing early maturing chickpea.


INTRODUCTION
Chickpea (Cicer arietinum L.) is a diploid annual legume with 2n = 16 chromosomes and a genome size of 738 Mb (Varshney et al., 2013). It is the world's second most important pulse crop after common bean with a total annual production of 13 million tons from an area of 13 million hectares (FAOSTAT, 2015). India, the largest producer and also the largest consumer of chickpeas in the world has 71% of global chickpea area.
Chickpea is a cool season crop mostly cultivated on residual soil moisture in the post-rainy season of the arid and semi-arid regions. Thus, the crop grows and matures on a progressively depleting soil moisture and experiences terminal drought (Kumar and Abbo, 2001). Terminal drought has become a major constraint in many chickpea growing areas. In addition, a large shift in chickpea area from cooler long-season environments to warmer short-season environments has increased the chances of exposure of crop to moisture and heat stresses at the reproductive stage causing severe yield losses . Early maturity has been identified as an important trait for increasing and stabilizing chickpea productivity by avoiding end of season drought (Subbarao et al., 1995;Kumar and Abbo, 2001) and frost (Warkentin et al., 2003) in short season environments. Significant impact of early maturing chickpea varieties in horizontal expansion of the crop in the semi-arid tropics has been reported (Than et al., 2007;Gaur et al., 2008).
Flowering time plays a key role in adaption and yield stabilization. It can be recorded with high precision and provides a good indication of succeeding phenological traits such as time of podding and maturity (Gaur et al., 2015). Large genotypic variations exist for flowering time in chickpea. Flowering time is a highly variable trait affected by various factors like soil moisture, photoperiod, temperature, altitude and latitude. The information available on genetics of flowering time in chickpea suggests that flowering time is controlled by one or a few major genes Kumar and van Rheenen, 2000;Anbessa et al., 2006;Hegde, 2010;Gaur et al., 2015;Gumber and Sarvjeet, 1996). Four flowering time genes have so far been identified in chickpea: efl-1 from ICCV 2 (Kumar and van Rheenen, 2000), efl-2 or ppd from ICC 5810 , efl-3 from BGD 132 (Hegde, 2010) and efl-4 from ICC 16641 (Gaur et al., 2015). Studies have shown that these flowering time genes are non-allelic (Hegde, 2010;Gaur et al., 2015).

Plant Material
Four early flowering lines ICCV 96029, ICC 5810, BGD 132 and ICC 16641 were crossed to a late flowering cultivar CDC Frontier. The early flowering lines chosen represent different sources of earliness genes based on the previous studies Kumar and van Rheenen, 2000;Hegde, 2010;Gaur et al., 2015). CDC Frontier is an improved kabuli cultivar with medium maturity developed at the University of Saskatchewan (Warkentin et al., 2005) and the genome sequence of this line was deciphered recently (Varshney et al., 2013). The F 1 seeds from all the crosses were selfed to develop F 2 mapping populations. The F 2 populations along with their parents and F 1 s were evaluated for segregation of flowering time during post rainy season of 2013-14 at International Crops Research Institute for the Semi-Arid Tropics (ICRISAT), Patancheru, Telangana, India. A total of 190 F 2 genotypes each of the crosses ICCV 96029 × CDC Frontier, ICC 5810 × CDC Frontier, BGD 132 × CDC Frontier and 146 F 2 s of the cross ICC 16641 × CDC Frontier were evaluated. F 3 progenies with 20 plants in each progeny row were evaluated for flowering time during post rainy season of 2014-15 with 164, 174, 182 and 102 progeny rows from each cross, respectively.

Phenotyping and Statistical Analysis
Flowering time was recorded as number of days from seeding to appearance of first fully opened flower on parents, F 1 s and F 2 populations. In F 3 , each progeny row was observed for flowering time at regular intervals and classified them as non-segregating and segregating types. Flowering time data was used to estimate the parameters of descriptive statistics and segregation analysis (Microsoft Excel, 2013). The expected values corresponding to the observed values for late and early flowering plants were calculated on the basis of the assumed Mendelian ratio. The deviations of these values were subjected to chi-square test to determine the goodness of fit.

Genomic DNA Extraction and Marker Genotyping
Young leaf tissues from 20 days old seedlings of parents and F 2 individuals were collected. Extraction of genomic DNA was carried out following high-throughput mini DNA extraction protocol as reported by Cuc et al. (2008). A total of 472 SSR markers reported earlier were used for parental polymorphism screening (including 146 CaM-series markers, Thudi et al., 2011;124 ICCM markers, Nayak et al., 2010;135 SSRs, Winter et al., 1999; 57 H-series markers, Lichtenzveig et al., 2005; and 10 NCPGR markers, Sethy et al., 2006;Gaur et al., 2011). SSR (CaM-series, ICCM-series, H-series, Winter-series and NCPGR markers) genotyping, PCR amplification, separation and visualization of amplified products were carried out by following the method described in previous studies (Nayak et al., 2010;Thudi et al., 2011).

Construction of genetic linkage map and QTL analysis
Genotypic data of all polymorphic markers from four mapping populations were compiled and linkage analysis was performed separately using JoinMap version 4.0 software (Van Ooijen, 2006) as described by Bohra et al. (2012). A consensus map was constructed based on data sets from four populations using JoinMap 4.0 software. QTL analysis of flowering time was carried out employing inclusive composite interval mapping (ICIM) using QTL-ICiMapping software version 4.0 (Wang et al., 2014). ICIM-Add mapping performs a stepwise regression to identify the most significant markers and marker-pair multiplications at 0.001 probability level and then scanning step of 1 cM. Later, a one-dimensional scanning or interval mapping was conducted to identify additive QTLs. The threshold levels to declare significance of a QTL was determined by performing 1,000 permutations by maintaining the chromosome-wise type I error rate of 0.05 (Churchill and Doerge, 1994). The LOD score peaks were used to estimate the most likely position of a QTL on the linkage map. The amount of variation explained by marker was determined using the coefficient of determination (R 2 ) value and expressed as percent phenotypic variance explained (PVE%). In this study, a QTL that explains more than 10% of total PVE was considered as major QTL.

Genetics of flowering Time genes
Flowering Time of Parental Lines, F 1 s and F 2 Populations The female parents ICCV 96029, ICC 5810, BGD 132 and ICC 16641 started flowering at 25, 28, 28 and 29 days after sowing, respectively. Whereas, the male parent CDC frontier flowered at 65 days (Supplementary Table 1

Individual Genetic Maps and Consensus Map
Among 472 SSRs tested, 100, 95, 90 and 93 were found polymorphic between the parents of the crosses ICCV 96029 × CDC Frontier, ICC 5810 × CDC Frontier, BGD 132 × CDC Frontier and ICC 16641 × CDC Frontier, respectively (Supplementary Table 4). Based on the distribution of markers on chickpea genome, a total of 76, 77, 68 and 68 polymorphic SSR markers were genotyped on the respective mapping population and then used for linkage map construction for each population separately. As a result, genetic linkage map for each cross and a consensus map were developed (http://cegresources. icrisat.org/cmap/sm/cp/mallikarjuna/) and the details are given below.

ICCV 96029 × CDC Frontier
A total of 75 SSR marker loci were mapped on 8 linkage groups (CaLGs) having a total map length of 248.76 cM and an average inter-marker distance of 3.32 cM (Table 2, Supplementary Figure 1). One marker (TA76s) remained unlinked to any of the linkage groups.

ICC 5810 × CDC Frontier
The genetic linkage map consisted of 75 SSR marker loci that are distributed across 8 linkage groups spanning 331.37 cM with an average marker density of one marker per 4.42 cM (Table 2, Supplementary Figure 2). Two markers (TA93 and TA76s) remained unlinked to any of the linkage groups.

BGD 132 × CDC Frontier
A total of 68 SSR marker loci were mapped on 8 linkage groups having a total map length of 311.10 cM and an average inter-marker distance of 4.58 cM ( Table 2,   *Null hypothesis of the test is that progeny segregate in the ratios tested. If the p-value (probability) is less than or equal to 0.05, then reject the null hypothesis. Otherwise one fails to reject the null hypothesis.
Supplementary Figure 3). No marker was found unlinked after linkage group assignment and ordering.

ICC 16641 × CDC Frontier
The intraspecific linkage map of this cross consisted of 67 SSR markers mapped onto 8 linkage groups spanning a total map length of 385.13 cM with an average marker density of 5.75 cM (Table 2, Supplementary Figure 4). Only one marker i.e., TA93 was unassigned to any of the linkage groups. Four genetic maps were integrated using JoinMap 4.0 to develop a consensus map that comprised of 8 linkage groups containing 109 markers with a total map length of 363.85 cM (Table 2, Figure 2). The map lengths of linkage groups in the consensus map were 32. 04,39.91,61.65,78.75,6.94,73.3,33.87 and 37.39 cM for CaLG01,CaLG02,CaLG03,CaLG04,CaLG05,CaLG06,CaLG07 and CaLG08 with 8,11,22,17,13,16,13 and 9 marker loci, respectively. The average density of the consensus map was 3.34 cM per marker.

QTLs for Flowering Time Genes
Ten significant QTLs were identified for flowering time in this study using QTL-ADD model of the ICIM software. The details of QTLs identified in each cross are presented below:

ICCV 96029 × CDC Frontier
A major QTL "Qefl1-2" was identified for flowering time on CaLG04 in the marker interval GAA47-ICCM192a with a peak LOD value of 5.66 and PVE of 11.75% (Table 3,  Supplementary Table 6, and Supplementary Figures 1, 5). In addition, a minor QTL "Qefl1-1" was also identified on CaLG03 in the marker interval CaM1122-TR13 with a peak LOD value of 3.45 and PVE of 5.66%. Both the QTLs showed negative additive effect indicating that allele for early flowering at this locus was contributed by ICCV 96029.

BGD 132 × CDC Frontier
One major and two minor QTLs were detected for the flowering time ( Table 3, Supplementary Table 6 and Supplementary Figures 3, 7). The major QTL "Qefl3-3" was located on CaLG08 in the marker interval TA127-H1D24. This was a highly significant QTL with a LOD peak value of 44.38 and PVE of 64.95%. The minor QTLs Qefl3-1 (LOD = 5.24; PVE = 4.39%) and Qefl3-2 (LOD = 4.21; PVE = 4.04%) were also detected on CaLG03 defined by marker intervals CaM1515-TR13 and TA142-TA64, respectively. In this cross also all the QTLs showed negative additive effect indicating that allele for early flowering at this locus was contributed by BGD 132.

ICC 16641 × CDC Frontier
A single major QTL (Qefl4-1) for flowering time was identified on CaLG06 flanked by markers TA14 and TR44 (Table 3, Figure 4B, Supplementary Table 6, and Supplementary Figure 4). This QTL contributed a high PVE of 88.19% at LOD value of 55.60. This QTL also showed a negative additive effect suggesting that allele for early flowering at this locus was contributed by ICC 16641. In this study, if more than one QTL of different cross share one or two flanking markers in common, it was considered as only one genomic region. The sequences of the markers flanking the QTL regions are provided (Supplementary Table 5).

DISCUSSION
Flowering time is an important trait for adaption of chickpea particularly in the semi-arid environments (Kumar and Abbo, 2001;Gaur et al., 2015). Information on genetic and molecular basis of flowering behavior would be useful for the breeding programs focusing on development of early maturing varieties. So far, four genes (efl-1, ppd/efl-2, efl-3 and efl-4) controlling flowering time in ICCV 96029 (Kumar and van Rheenen, 2000), ICC 5810 Hegde, 2010), BGD 132 (Hegde, 2010) and ICC 16641 (Gaur et al., 2015) have been reported in chickpea. When these early flowering lines were crossed with a late flowering cultivar (CDC Frontier), F 1 s were late to flower as the gene for delayed flowering is known to be dominant to early flowering in chickpea Kumar and van Rheenen, 2000;Hegde, 2010;Gaur et al., 2015). Segregation analysis (  Therefore, the present study confirmed the single recessive gene hypothesis for flowering time in ICCV 96029 (Kumar and van Rheenen, 2000), BGD 132 (Hegde, 2010) and ICC 16641 (Gaur et al., 2015) and the digenic control in ICC 5810 (Hegde, 2010;Gaur et al., 2015). This implies that the early flowering trait can be easily incorporated into the desired genetic backgrounds.
In the present study though sufficient number of SSRs (472) that represent most of the chickpea genome were used, a low polymorphism level (21.40, 20.13, 19.07 and 19.70%) was observed compared to previous studies (Tar Kottapalli et al., 2009). The goodness of fit of the observed segregation ratio to the expected ratio demonstrated that the majority of the SSRs did not significantly deviate from the expected 1:2:1 ratio (P ≥ 0.05). Since F 2 populations were used, the intraspecific maps developed in this study represent the coarse genetic maps spanning 248.58, 331.37, 311.10 and 385.13 cM in ICCV 96029 × CDC Frontier, ICC 5810 × CDC Frontier, BGD 132 × CDC Frontier and ICC 16641 × CDC Frontier, respectively. These results indicate that the intraspecific maps obtained are less dense compared to earlier maps (Winter et al., 2000;Nayak et al., 2010;Thudi et al., 2011;Varshney et al., 2014). Also, a varying levels of marker distributions were observed with dense sub-clusters of marker loci either in the central region or at distal ends of most of the linkage groups in all the maps. It may reflect the low level of recombination in centromeric and subtelomeric genomic regions (Tanksley et al., 1992) and such apparent clustering of markers on the linkage groups was also observed in the previous studies (Tanksley et al., 1992;Winter et al., 1999;Millan et al., 2010;Nayak et al., 2010). When these genetic maps were compared with earlier maps, the SSRs that are common between the current maps and the previous maps (Millan et al., 2010;Thudi et al., 2011;Jaganathan et al., 2014;Varshney et al., 2014) were placed on the same linkage groups which encourages the possibility of integration of different maps through common markers. However, the order of marker loci on intra-specific maps differed in several instances.
Construction of a consensus map based on synteny between the several linkage maps will make it possible to expand the chickpea genetic map and increase the marker density. In the present study, the data sets from four populations were joined to develop the consensus map. While comparing four intraspecific genetic maps, 28 loci were found common between four maps. These were considered as anchor markers and used for merging the genetic maps for construction of the consensus genetic map. The consensus map developed contained 109 markers that covered 363.85 cM of map length, which is less dense compared to the consensus maps of Millan et al. (2010) and Varshney et al. (2014). A detailed comparison between individual component maps and consensus map reflects a general coincidence. Although differences in marker order exist, linkage groups are generally conserved. The sub clusters and gaps were also observed in most of the LGs either at central or in distal FIGURE 3 | Mapping of major flowering time gene "efl-3" on CaLG08 of the cross BGD 132 × CDC Frontier. (A) Mapping of major flowering time gene "efl-3" on CaLG08 based on F 3 segregating data (B) Mapping of major QTL for flowering time "Qefl3-3" on CaLG08.
regions. Since the consensus map is low in marker density, saturating these regions with more markers will help in any map based cloning of agronomically important genes.
In the present study, 10 genomic regions were identified for flowering time including seven QTLs having major effects with PVE more than 10%. A major QTL "Qefl1-2" (PVE = 11.75%) for flowering time was detected in ICCV 96029 × CDC Frontier on CaLG04 (flanked by GAA47-ICCM192a). Mendelian inheritance revealed that flowering time was governed by a single major gene. The identified QTL on CaLG04 could be same as the chromosomal region reported by Daba et al. (2016) who mapped four QTLs for days to flowering on LG4 in the same cross. Previously, Cobos et al. (2007) reported a major QTL for days to 50% flowering (QTL DF1 ; 20% PVE) on LG4. This QTL had a common marker (i.e., GAA47) with the QTL reported in this study. Therefore, they may refer to the same QTL. In the present study, a minor QTL "Qefl1-1" (PVE = 6.90%) was also identified on CaLG03. Recently, Daba et al. (2016) also mapped a minor QTL on LG3 in the same cross. Therefore, these QTLs may be representing the same genomic region. Cho et al. (2002) and Jamalabadi et al. (2013) also reported a QTL for days to flowering on LG3 using a RIL population from a cross involving ICCV 2 as one of the parents. The line ICCV 2 is an indirect source of earliness (efl-1) in our study. However, lack of common markers does not allow a definitive conclusion that these two QTLs represent the same locus. Daba et al. (2016) identified additional major QTLs for days to flowering on LG5 and LG8 that are consistent across years and sites. However, no such QTLs were identified in the present study since we evaluated the mapping population at only one location. Based on these findings it is apparent that several unknown factors confer time to flowering in chickpea even though segregation for a major flowering gene was observed. Similar findings were also reported earlier (Cho et al., 2002). In ICC 5810 × CDC Frontier, major QTLs on CaLG01 (Qefl2-1, PVE = 20.28%), CaLG03 (Qefl2-2, PVE = 24.95%), CaLG04 (Qefl2-3, PVE = 10.53%) and CaLG08 (Qefl2-4, PVE = 25.73%) were identified. Genetic studies revealed that two major genes with complementary gene action controlling flowering time in this cross. Earlier, Cho et al. (2002) detected a single QTL for flowering time on LG3 between the markers TS57 and TA127. Recently, Jamalabadi et al. (2013) also identified a QTL for flowering time on LG3 closely linked to the marker TA117. However, these markers were not mapped in the present study and hence the exact chromosomal location could not be compared. Whereas, Cobos et al. (2009) andAryamanesh et al. (2010) mapped a QTL for flowering time on LG3 closely linked to marker TA142 which was also detected in this study. Therefore, these QTLs could belong to the same set of genes. Another QTL for flowering time was identified by Cobos et al. (2007) on LG4 (explaining 20% PVE) closely linked to the marker GAA47. A QTL on CaLG4 (Qefl2-3, PVE = 10.53%) was detected having GAA47 as flanking marker in the present study. Therefore, these QTLs could be same in both the findings. In all these studies, different parental lines were used. However, Lichtenzveig et al. (2006) in the cross involving Hadas and ICC 5810 reported three QTLs on LG1, LG2 and LG8 for flowering time. Recently, Rehman et al. (2011) reported four QTLs for flowering time on LG1, LG3, LG4 and LG8. None of these QTLs were found similar to the QTLs detected in the present study. In our study, however, LG2 was not associated with any effect on flowering time. One possible explanation for this could be the absence of common markers in our map due to non-availability of more number of polymorphic markers for linkage analysis. Further studies are needed to confirm which two major genes out of four QTLs detected in this study are responsible for flowering time in ICC 5810.
In BGD 132 × CDC Frontier, a major QTL "Qefl3-3" (flanked by markers TA127 and H1D24) was detected for flowering time on CaLG08 with 64.95% PVE (Figure 3B). This is the first report of mapping major flowering time gene "efl-3" in BGD 132. Linkage analysis based on F 3 segregating data resulted in mapping of flowering time locus "efl-3" on CaLG08 in this cross ( Figure 3A). Previously, Cho et al. (2002) reported a QTL for flowering time on LG3 flanked by markers TS57 and TA127. However, LG3 of Cho et al. (2002) is equivalent to CaLG8 in this study based on the common markers of the current map and genetic maps of Tar'an et al. (2007) and Varshney et al. (2014). It appears that both the QTLs could be the same. Two additional minor QTLs i.e., "Qefl3-1" and "Qefl3-2" were also detected on CaLG03. Hence, CaLG08 appears to be a strong candidate linkage group having major QTL controlling flowering time gene "efl-3." In ICC 16641 × CDC Frontier, a single putative QTL (Qefl4-1; PVE = 88.19%) for flowering time was detected on CaLG06 between markers TA14 and TR44 ( Figure 4B). This novel QTL is unique for the flowering time gene "efl-4" and is reported for the first time in this study. Mendelian inheritance also revealed monogenic inheritance of flowering time in this cross. This was further confirmed by linkage analysis and mapping of major flowering time gene "efl-4" on CaLG06 of the chickpea genetic map between the markers TA14 and TR44 ( Figure 4A). Therefore, this genomic region can be targeted for developing early maturing chickpea varieties through Marker Assisted Breeding (MAB).

CONCLUSIONS
The present study revealed major gene inheritance of flowering time genes under short season environment typical to semi-arid tropics. This simple inheritance of early flowering trait can be easily transferred into desired genetic backgrounds. The SSRs were used to construct separate intraspecific linkage maps of four F 2 populations. These linkage maps were combined to construct a consensus map of chickpea with 109 marker loci (363.85 cM). QTL analysis revealed altogether 10 genomic regions for flowering time including seven major QTLs distributed across CaLG01, CaLG03, CaLG04, CaLG06 and CaLG08 of chickpea genetic map. It is also the first report on mapping of flowering time genes "efl-3" and "efl-4" in chickpea. These genomic regions provide strong basis for further investigation on fine mapping and validation of the identified QTLs which will help in developing early-maturing chickpea varieties under short season environments.