Quantitative Trait Loci for Cold Tolerance in Chickpea

Fall-sown chickpea (Cicer arietinum L.) yields are often double those of spring-sown chickpea in regions with Mediterranean climates that have mild winters. However, winter kill can limit the productivity of fall-sown chickpea. Developing cold-tolerant chickpea would allow the expansion of the current geographic range where chickpea is grown and also improve productivity. The objective of this study was to identify the quantitative trait loci (QTL) associated with cold tolerance in chickpea. An interspecific recombinant inbred line population of 129 lines derived from a cross between ICC 4958, a cold-sensitive desi type (C. arietinum), and PI 489777, a coldtolerant wild relative (C. reticulatum Ladiz), was used in this study. The population was phenotyped for cold tolerance in the field over four field seasons (September 2011–March 2015) and under controlled conditions two times. The population was genotyped using genotypingby-sequencing, and an interspecific genetic linkage map consisting of 747 single nucleotide polymorphism (SNP) markers, spanning a distance of 393.7 cM, was developed. Three significant QTL were found on linkage groups (LGs) 1B, 3, and 8. The QTL on LGs 3 and 8 were consistently detected in six environments with logarithm of odds score ranges of 5.16 to 15.11 and 5.68 to 23.96, respectively. The QTL CT Ca-3.1 explained 7.15 to 34.6% of the phenotypic variance in all environments, whereas QTL CT Ca-8.1 explained 11.5 to 48.4%. The QTLassociated SNP markers may become useful for breeding with further fine mapping for increasing cold tolerance in domestic chickpea. D. Mugabe, Dep. of Crops and Soil Sciences, Washington State Univ., Pullman, WA 99164; C.J. Coyne and E. Landry, USDA-ARS Plant Germplasm Introduction and Testing, Washington State Univ., Pullman, WA 99164; J. Piaskowski, College of Agriculture and Life Sciences, Univ. of Idaho, Moscow, ID 83844; P. Zheng, Y. Ma, and D. Main, Dep. of Horticulture, Washington State Univ., Pullman, WA 99164; R. McGee and G. Vandemark, USDA-ARS Grain Legume Genetics and Physiology, Pullman, WA 99164; H. Zhang, Dep. of Soil and Crop Sciences, Texas A & M Univ., College Station, TX 77843; S. Abbo, Faculty of Agriculture, Hebrew Univ. of Jerusalem, Rehovot, Israel. Received 17 Aug. 2018. Accepted 3 Nov. 2018. *Corresponding author (clarice.coyne@ars.usda. gov). Assigned to Associate Editor Manjit Kang. Abbreviations: CF, Central Ferry farm; CIM, composite interval mapping; FT1, Freezing Test 1; FT2, Freezing Test 2; GBS, genotypingby-sequencing; LG, linkage group; LOD, logarithm of odds; QTL, quantitative trait locus/loci; RIL, recombinant inbred line; SNP, single nucleotide polymorphism. Published in Crop Sci. 59:573–582 (2019). doi: 10.2135/cropsci2018.08.0504 © Crop Science Society of America | 5585 Guilford Rd., Madison, WI 53711 USA This is an open access article distributed under the CC BY license (https:// creativecommons.org/licenses/by/4.0/). Published February 14, 2019

diversification.As a spring-sown crop, the short growing season of chickpea limits its grain yield, leaves very little crop residues to combat soil erosion, and hardly contributes to soil organic matter (Kumar and Abbo, 2001).The development of a fall-sown, winter-hardy chickpea could improve yields by extending the growing season, help the crop escape heat stress during flowering and late-season drought, and provide additional protection from cool temperatures during seed set in the early spring (Singh et al., 1997).Chickpea lacks winter hardiness because of genetic bottlenecks and subsequent loss of critical alleles during its evolution under domestication (Abbo et al., 2003).However, alleles for winter hardiness and a vernalization requirement still exist and are prevalent in chickpea's wild progenitor, C. reticulatum Ladiz (Abbo et al., 2002).
Chickpea wild relatives have been recommended as sources of cold tolerance for decades (van der Maesen and Pundir, 1984;Singh et al., 1998;Toker, 2005).In a study that compared wild with domesticated chickpea, it was reported that wild species had significantly higher levels of cold tolerance than the best cold-tolerant cultivars (Toker, 2005).Singh et al. (1998) concluded that C. reticulatum was one of the most cold-tolerant annual Cicer species.The use of crop wild relatives, such as C. reticulatum, carries significant agronomic challenges, such as reticulated seed, prostrate growth habit, and pod dehiscence (Ladizinsky, 1975;Ladizinsky and Adler, 1976).However, C. reticulatum holds the advantage of a vernalization requirement that contributes to winter hardiness by delaying flowering (Abbo et al., 2002).Using interspecific recombinant inbred line (RIL) populations, vernalization in chickpea was recently mapped to Linkage Group (LG) 3 (Samineni et al., 2016;van-Oss et al., 2018).
Understanding the genetic basis of cold tolerance in chickpea is important to be able to introgress novel alleles efficiently into elite germplasm.Thus, the objectives of the current study were to phenotype and genotype a RIL population and identify QTL associated with cold tolerance in chickpea.The identified QTL-associated single nucleotide polymorphism (SNP) markers have potential for use in introgressing improved cold tolerance into domesticated chickpea.

Plant Material
The chickpea population used in this study was 127 F 6 -derived F 10 RILs (CRIL2, C = chickpea, 2 = number of populations in the USDA breeding program) from the wide cross of C. arietinum ICC4958 ´ PI 489777 C. reticulatum.PI 489777 C. reticulatum is the wild progenitor of domestic chickpea that has a prostrate growth habit and small (15.4 g, 100-seed weight), dark brown seeds covered by tiny tubercules (Ladizinsky and Adler, 1976).Cultivar ICC4958 is a semi-spreading desi type with larger (34.0 g, 100 seed weight), brown seeds (Saxena et al., 1993).PI 489777 is the cold tolerance source and ICC4958 is the cold-susceptible parent.The cross was made by the USDA in 1992 (X92C031) to study Fusarium wilt resistance (caused by Fusarium oxysporum Schlechtend.),and RILs were created via single-seed descent in a greenhouse (Brim 1966;Snape et al., 1985;Ratnaparkhe et al., 1998;Tekeoglu et al., 2000).

Controlled Conditions
Controlled environment experiments (Freezing Test 1 [FT1] and Freezing Test 2 [FT 2]) were conducted using 104 RILs because of lack of sufficient seed, and the two parents of CRIL2.The experimental design for both experiments was a randomized complete block design repeated twice, where each variants were identified using the genotypes program in Stacks (Catchen et al., 2013).A total of 747 SNP markers were used to construct a linkage map (Supplemental Table S5).A maximum of 0.2 was set for the allowable missing portion of missing data for SNPs.A linkage map was constructed using Joinmap 4.0 (Van Ooijen, 2006), and minimum logarithm of odds (LOD) for group selection was five, with recombination frequency set between 0.25 and 0.05 cM.Distances were calculated using the Haldane's mapping function (Haldane, 1919).The Monte Carlo maximum likelihood algorithm (Geyer and Thompson, 1992) was used for ordering the markers.The segregation ratio at each marker locus was statistically analyzed against the expected Mendelian segregation ratios of 1:1 using c 2 tests.

Data Analysis
The data were analyzed separately for each response variable (percentage leaf damage or percentage surviving stand count).Genotype and replication (blocking factor) were treated as random effects, and year was treated as a fixed effect.For freeze tests, genotype was treated as a random effect, and replication and experiment were treated as fixed effects, as each had only two levels.
The statistical model used was where ijk Y is the response variable, b 0 is the overall mean, and b 1 is the effect of year.The random variables u 1 , u 2 , and e ijk (the effects from replication, the genotype effects, and the error term, respectively) were distributed as follows: The data (for leaf damage and stand count) being binomial traits, a generalized linear model was used for each variable using the binomial family and logit link function: The data were also analyzed with a general linear model using the square root transformation, log transformation, and arcsine of the square root transformation on the response data.The freeze test data were analyzed using a general linear model.Residuals from all models were examined for homoscedasticity and normality.The experiment was analyzed in R version 3.3.4with Asreml (Butler et al., 2007;R Development Core Team, 2015).Heritability estimates were generated for each environment using the Sommer package in RStudio 2015 (Covarrubias-Pazaran, 2018).Broad-sense heritability was calculated as replication entry had four plants.Each freezing test consisted of 104 RILs, eight pairs of the parents and check lines.Based on 2014 and 2015 field data, the five most cold-tolerant RILs (CRIL2-111, -96, -122, -53, and -5) and the five most susceptible RILs (CRIL2-1, -11, -85, -29, and -47) were selected as checks.The freezing protocol used in the freezing chamber was based on the procedure for winter wheat (Zhu et al., 2014).Seeds were scarified before each experiment using a multitool (Dremel Multipro, Ball Bearing Model 5) and germinated on moist paper at 20°C for 4 d, then planted in a soilless mix Sunshine #5 (Sun Gro Horticulture) in 6.35-cm ´ 6.35-cm ´ 5.66-cm trays.Four plants per replication were grown for 3 wk in the greenhouse maintained at 22 to 25°C/16 to 19°C day/ night temperatures.Slow-release fertilizer Osmocote (15-9-12 N-P-K) with micronutrients (The Scotts Company) was applied 2 wk after planting.Plants were then moved to a growth chamber for 5 wk.Temperatures in the growth chamber were recorded using a HOBO data logger (Onset Manufacturing Company).Acclimation started at 1 wk with 10°C and 16 h of light, 2 wk with 5°C and 12 h of light, and 2 wk with 3°C and 10 h of light.The last step was freezing that started at −3°C with no light for 16 h, then −5°C with no light for 2 h.Temperature was dropped further to −7°C with no light for 2 h during FT 1 because of a freezing chamber technical issue.From the freezing growth chamber, the plants were moved to the acclimation growth chamber at 4°C for 24 h, then returned to initial greenhouse conditions, and the plants were rated for cold damage after 7 d.Evaluation of cold damage was based on leaf damage that was scored using a 1-to-9 visual scale developed by Fiebelkorn (2013) (Table 1, Supplemental Tables S4a and S4b).

Linkage Map Construction
DNA was extracted from leaf samples collected from 124 3-wkold, F 6 -derived F 10 RILs using the Qiagen DNeasy 96 Plant Kit.Genotyping-by-sequencing (GBS) libraries were prepared using the two-enzyme method described by Poland et al. (2012) that uses methylation-sensitive restriction enzymes (PstI/MspI).Libraries were sequenced using Illumina Hiseq 2500.The Illumina DNA-sequencing raw reads were demultiplexed and cleaned using the process-radtags program in Stacks (Version 1.47;Catchen et al., 2013), and high-quality reads were aligned on reference genome developed from the kabuli cultivar 'CDC Frontier' (Varshney et al., 2013), as per Thudi et al. (2016), using the Burrows Wheeler Alignment tool (Version 0.5.9) with default alignment parameters (Li and Durbin, 2009).Sequence where V G is the genotype variance, V GE is the genotype ´ environment variance, V E is the environmental variance, and N env is the number of environments.

QTL Detection with Composite Interval Mapping
A genetic linkage map made up of 747 SNPs and cold tolerance phenotypic data obtained from each of the six environments was used for QTL analysis using "R/qtl" in R software version 3.3.4by employing the composite interval mapping (CIM) method (Broman et al., 2003(Broman et al., , 2010;;Arends et al., 2010).The statistical significance of the QTL was assessed by using Churchill and Doerge's permutation test with 1000 replications and a significance level of P £ 0.05 (Churchill and Doerge, 1994).The additive effects and R 2 of the detected QTL were estimated using the "fitqtl" function of R version 3.3.4.The 95% Bayes interval was used to obtain interval estimates of QTL locations.

RESULTS
As expected, the two parents consistently differed in their reactions to the winter temperatures, and transgressive segregants with increased cold tolerance were identified in most environments (Fig. 1, Supplemental Tables S2 and S3).The 2011-2012 field study had the mildest low temperatures (Supplemental Fig. S1A), and no plant death was observed.The subsequent 3 yr reached lower temperatures (−12.5, −5.9, and −17°C, respectively; Supplemental Fig. S1B-S1D), with accompanying plant death of the most susceptible lines and the C. arietinum parent ICC 4958.
Given the residual pattern, field leaf damage (2011)(2012)(2012)(2013), and stand count were best modeled with the binomial distribution and the logit link function.Based on the residual pattern, stand count (2013)(2014)(2014)(2015) was modeled after an arcsine transformation using a standard linear model with normal distribution.The variance components indicated that genotypes explained a high percentage of the total variance in the 4 yr of field studies but a low percentage of the variance in the freezing tests (Table 2).Additionally, there was a relatively high error term in the freeze tests.
The predicted values for leaf damage and stand count were moderately correlated with each other (r = 0.62, r = 0.64 [Pearson and Spearman rank correlations, respectively]).Replication consistently had no effect, as either a fixed effect or a random effect.Correlations between leaf damage years and stand count years were higher (r = 0.75 and 0.85, respectively) using untransformed data (Supplemental Table S1).Correlations between the 4 yr of field cold tolerance studies were higher than between the two freeze tests (FT2, Table 3).Heritability estimates were as follows: CF 2012 and 2013 had an H 2 of 0.82 (SE = 0.03), and CF 2014 and 2015 had an H 2 of 0.96 (SE = 0.01) for the field environments, and controlled environments had an H 2 of 0.67 (SE = 0.07).

Linkage Map
An interspecific genetic linkage map containing 747 SNPs mapped on eight LGs (Supplemental Table S5, Supplemental Fig. S2) was constructed.The map spanned a total distance of 393.7 cM at an average density on 1.8 SNPs cM −1 .Nearly half of the loci (48%) exhibited segregation distortion from the expected ratio of 1:1 for RILs.The smallest LG in size was LG 7, which spanned 24.8 cM and contained 58 markers.Linkage Group 4 was the largest, with 203 markers spanning 90.3 cM.The densest LG, LG 5 (A), had an average marker density of 2.8 SNPs cM −1 .

QTL detection using Composite Interval Mapping
Using CIM QTL analysis for cold tolerance, three QTL were identified, of which two were large-effect QTL that were significant in all six environments, whereas one QTL was only detected in one environment (Table 4, Supplemental Table S6).The minimum LOD score set by Churchill and Doerge's permutation test (Churchill and Doerge, 1994) was 3.1.All the detected QTL for cold tolerance had C. reticulatum (PI 489777) as the contributing parental allele.The QTL were named following Hamon et al. (2011).The first QTL, CT Ca-3.1, was detected on LG 3 at 43.8 cM in all the six environments.CT Ca-3.1 had a confidence interval between 0.5 and 13.4 cM.The minimum phenotypic variation for cold tolerance trait explained by this QTL was 7.15%, with a LOD score of 6.89 in the 2015 CF field environment, whereas the maximum was 34.57%, with a LOD score of 15.11 in the 2013 CF field environment.The second QTL, CT Ca-8.1, was detected on LG 8 in the four CF field sites and two controlled environments.The confidence interval for CT Ca-8.1 was between 0 and 9.04 cM.The minimum phenotypic variation for cold tolerance trait explained by CT Ca-8.1 was 11.47%, with a LOD score of 5.68 in FT2, whereas the maximum was 48.41%, with a LOD score of 23.96 in the 2014 field environment.The third QTL, CT Ca-1.1, was detected on LG 1(B) in the 2015 CF field environment.The phenotypic variation for cold tolerance trait explained by this QTL was low (2.6%), with a LOD score of 3.14, just above the threshold to declare a QTL.
The minimum LOD score was set at 3.4 for plant height and seed volume, and significant QTL were identified (Table 4).Two QTL were identified for plant height: Ht Ca-4 was found on LG 4, with a LOD score of 6.5, whereas Ht Ca-8 was found on LG 8, also with a LOD score of 6.5.Ht Ca-4 explained 20.21% of the phenotypic variance, whereas Ht Ca-8 explained 19.97% of the phenotypic variance.For seed ellipsoid volume, two QTL were identified.The first QTL, SEV Ca-1(A), was found on LG 1(A), with a LOD score of 7.4, explaining 15.93% of phenotypic variance.The second QTL, SEV Ca-4, was found on LG 4, with a LOD score of 11.8, explaining 29.41% of phenotypic variance.tolerance was explained by CT Ca-3.1 and CT Ca-8.1, the QTL on LGs 3 and 8, respectively.As expected, all largeeffect QTL that were associated with cold tolerance were derived from C. reticulatum.Similar results were reported in pea, where a QTL (gene region Hr) explained up to 43% of the variation for cold tolerance (Lejeune-Hénaut et al., 2008;Klein et al., 2014).It is also similar for barrel medic, lentil, and faba bean, where the one consistent QTL explained up to 39, 28, and 23% of the variation, respectively (Kahraman et al., 2004;Avia et al., 2013;Sallam et al., 2017).This is not surprising, given the moderate to very high broad-sense heritability (0.67-0.96) for cold tolerance in the CRIL2 population.Four consistent cold

DISCUSSION
The two most significant QTL, identified using CIM, for cold tolerance in the CRIL2 chickpea population explained high percentages of the phenotypic variance that could be used in breeding winter-hardy chickpea.Up to 34 and 48% of the phenotypic variance for cold Table 2. Variance components of the field study and the controlled freeze tests.
Although the Mallikarjuna et al. (2017) chickpea consensus linkage map (363.8 cM) is similar to this report's linkage map (393.7 cM, Supplemental Table S5), there are gaps (Supplemental Fig. S2) compared with high-density chickpea maps (Deokar et al., 2018).We hypothesize that the two enzymes, PstI/MspI, of the Poland et al. (2012) GBS were not the best combination for chickpea, since only 747 SNPs were identified.Currently, we are testing the one enzyme ApeKI GBS of Elshire et al. (2011) for chickpea.The SNPs in CRIL2 were successfully identified using CDC Frontier as the reference genome, which has emerged as a very useful and highly cited chickpea reference genome (Thudi et al., 2016).Almost half of the markers exhibited segregation distortion and were not used in the linkage map.Hybrid progeny populations derived from interspecific crosses between domesticated cultivars and wild accessions of their immediate wild progenitors usually show higher polymorphism levels as compared with intraspecific crosses between domesticated parents.Such polymorphism is fundamental for our ability to produce accurate and dense genetic maps.However, a prominent weakness of such experimental populations is the distorted, non-Mendelian allelic segregation, which is an almost universal phenomenon of genetic analyses based on interspecific cross progeny, as evident from numerous Table 3. Pearson phenotypic correlation coefficients between the cold tolerance ratings determined in the four field studies (percentage damage) and two freeze tests (FT, plant damage scores) for cold tolerance of the ICC4958 ´ PI 489777 CRIL2 population (Supplemental Tables S1, S4a, and S4b).reports that need no reiteration here.Indeed, for many decades, geneticists have been lamenting the occurrence of unequal segregation of alleles in such crosses (Zamir and Tadmor, 1986) and likewise the differential survival of zygotes after interspecific crosses (Gadish and Zamir, 1987).For this reason, the high frequency of distorted segregation observed in the present work is no exception.
Although CT Ca-3.1 and CT Ca-8.1 cold tolerance QTL were detected with the freeze test (FT1 and FT2) data, the high experimental errors and low correlation values indicated that the testing procedure for chickpea would need to be improved.Assessing leaf damage at different temperature drop points and rates is important for improving the procedure in future studies.Assessment of additional traits such as plant roots, tissues, turgidity, and discoloration after each freeze test may capture the full severity of plant cold damage not detected using just the visual scores of Table 1 (Landry and Hu, 2018).The ability to complete a cold tolerance test every 8 wk throughout the year as opposed to once per year in a field study will allow breeders to select more frequently and cycle generations faster, thus decreasing the amount of time required for variety development.
The QTL CT Ca-3.1 is a good candidate region for the vernalization gene previously mapped in the CRIL2 population and in one other population where a C. reticulatum accession was one of the parents (Samineni et al., 2016;van-Oss et al., 2018).The same QTL was mapped in this study.Vernalization is the acceleration of flowering after a plant's exposure to cold temperature (Sung and Amasino, 2004).Flowers are known to be delicate structures on a plant that are sensitive to various stresses (Kim et al., 2009), and vernalization protects a plant from cold stress by stopping it from flowering before winter and allowing it to flower when temperatures are optimal (Kim et al., 2009).The vernalization requirement has been well described in C. reticulatum (Abbo et al., 2002;Sharma and Upadhyaya, 2015;van-Oss et al., 2015van-Oss et al., , 2018)).We mapped the vernalization gene (MtVRN2), a repressor of the flowering locus T (FT) gene homolog from Medicago truncatula, to the CTCa-3.1 confidence interval using the CDC Frontier chickpea reference genome ( Jaudal et al., 2016).This is a start at candidate gene identification; fine mapping will be necessary to elevate its status.A draft genome was published for both CRIL2 parents, which should facilitate this endeavor ( Jain et al., 2013;Gupta et al., 2016).In a more distantly related legume, Lupinus angustifolius L., a gene (LanFTc1) also interacts with flowering locus T (FT), controlling its vernalization requirement in lupin (Nelson et al., 2017) and thus rendering the M. truncatula gene especially interesting in chickpea.
The frost tolerance gene Hr in pea (early flowering 3, ELF3 ortholog) is also a flowering control gene of interest (Lejeune-Hénaut et al., 2008;Weller et al., 2012).Recently, Efl1, an ELF3 ortholog in chickpea, was mapped to the early flowering QTL on chickpea LG5, so it likely is not the vernalization gene in chickpea (Ridge et al., 2017).Other genes of interest are the CBF/DREB1 genes, which are key to cold acclimation in numerous plants, with evidence published for their role in pea and barrel medic cold tolerance (Thomashow, 2010;Tayeh et al., 2013).A Basic Local Alignment Search Tool (BLAST) search revealed a homolog on LG 4 of chickpea, but it was not found in this QTL study, possibly due to linkage map density of LG 4 (Supplemental Table S5, Supplemental Fig. S2).
Several QTL for chickpea plant height have been identified in other studies on chickpea (Gowda et al., 2011;Kujur et al., 2016).Although two studies reported QTL for plant height on the same LGs as in this study (Hgt Ca-4 and Hgt Ca-8), they were not in the same regions, ranging from 31.3 to 37.3 cM distal to Hgt Ca-4 and Hgt Ca-8, respectively (Gowda et al., 2011;Kujur et al., 2016).These two studies used C. arietinum intraspecific RIL populations that did not segregate for the prostrate growth habit of C. reticulatum.Therefore, this was not unexpected.However, two studies, using a similar wide cross, mapped a genegene region for prostrate growth habit to LG 3 (Cobos et al., 2009;Aryamanesh et al., 2010) that was not detected in this study.Further, Upadhyaya et al. (2017) identified eight SNPs on LGs 1, 3, 4 (3), 6, and 7 of Scaffold682 associated with plant growth habit using categories rather than height per se.The QTL on LG 4 does co-localize with Ht Ca-4.
For chickpea seed size, many studies have reported QTL by using 100-seed weight (Cobos et al., 2007;Gowda et al., 2011;Karami et al., 2015).While none of these studies reported QTL for seed size using seed volume, all three reported QTL for seed size on LG 4, where the SEV Ca-4 seed volume QTL resides.The closest published QTL of Gowda et al. (2011) was 3 cM from the SEV Ca-4 marker, well within the confidence interval to be considered potentially the same QTL.
The identification of these QTL is a start on efficiently moving positive alleles for cold tolerance into cultivated chickpea.The main QTL identified in this study have large confidence intervals (13.4 and 9 cM).Creating a higher density map and performing fine mapping may help reduce the confidence interval and make SNP discovery more precise (Gao and Zhu., 2013;Jaganathan et al., 2015).This can increase the efficiency of moving the cold tolerance trait from C. reticulatum into domesticated chickpea.

Fig. 1 .
Fig. 1.Frequency histogram of the cold tolerance ratings of the CRIL2 recombinant inbred line (RIL) population across 4 yr in the field, measured by percentage leaf damage or percentage stand count, and the two controlled freeze tests measured with a 1-to-9 ordinal scale, with 9 being dead plants.P1 is PI 489777, the cold-tolerant parent, whereas P2 is ICC 4958, the cold-susceptible parent.
2011 on the USDA Central Ferry (CF) farm; 28 Sept. 2012 on the CF farm and Washington State University Whitlow and Spillman farms; 25 Sept. 2013 on the CF farm and Spillman farm; and 1 Oct. 2014 on the CF farm.

Table 1 .
Fiebelkorn (2013)sual scores assigned to CRIL2 plants in the controlled environment experiments of Freeze Test 1 and 2 as developed byFiebelkorn (2013).