Genome-Wide Association Mapping Reveals Multiple QTLs Governing Tolerance Response for Seedling Stage Chilling Stress in Indica Rice

Rice crop is sensitive to cold stress at seedling stage. A panel of population representing 304 shortlisted germplasm lines was studied for seedling stage chilling tolerance in indica rice. Six phenotypic classes were exposed to six low temperature stress regimes under control phenotyping facility to investigate response pattern. A panel of 66 genotypes representing all phenotypic classes was used for ensuring genetic diversity, population structure and association mapping for the trait using 58 simple sequence repeat (SSR) and 2 direct trait linked markers. A moderate level of genetic diversity was detected in the panel population for the trait. Deviation of Hardy-Weinberg's expectation was detected in the studied population using Wright's F statistic. The panel showed 30% variation among population and 70% among individuals. The entire population was categorized into three sub-populations through STRUCTURE analysis. This revealed tolerance for the trait had a common primary ancestor for each sub-population with few admix individuals. The panel population showed the presence of many QTLs for cold stress tolerance in the individuals representing like genome-wide expression of the trait. Nineteen SSR markers were significantly associated at chilling stress of 8°C to 4°C for 7–21 days duration. Thus, the primers linked to the seedling stage cold tolerance QTLs namely qCTS9, qCTS-2, qCTS6.1, qSCT2, qSCT11, qSCT1a, qCTS-3.1, qCTS11.1, qCTS12.1, qCTS-1b, and CTB2 need to be pyramided for development of strongly chilling tolerant variety.


INTRODUCTION
Rice is life for millions of people in the globe. Globally, rice the most important food crop for 3.5 billion populations with a total production of 715 million tons (FAO, 2013). Rice cultivation covers around 160 million hectares on the earth. In India, around one-fourth (43 million hectares) of the world's rice area is cultivated for rice with a total production of 104.8 million tons (DAC, 2015). Globally, around 15 million hectares of rice fields are affected by cold weather spreading in 24 different countries. Currently, there is a ceiling of rice yield potential which is not sufficient to feed the ever growing population. Therefore, it is essential to increase the rice production even from abiotic stressful rice-growing ecosystems to enable the increasing food demand to commensurate with increasing population. Around 4 million hectares of rice areas comprising boro and part of dry season rice of India are highly affected by seedling-stage cold causing delay in growth of the plant which subsequently further aggravated as it coincides with high temperature during flowering stage . The ideal temperature for rice germination and seedling growth is 25-35 • C and the temperature below 15 • C normally affects various processes (Nakagahra et al., 1997). Rice plants at early vegetative stage are sensitive to cold stress resulting in poor germination, slow early growth, yellowing and withering of plant, reduced tillering, and stunted growth which finally decrease in overall yield (Zhang et al., 2005;Andaya and Tai, 2006;Lou et al., 2007;Suh et al., 2010;Pradhan et al., 2015). The extent of cold and high temperature stress is increasing due to global climate change. For instance, in Asia sudden low temperature and high temperature stress covering areas have been widened (Pradhan et al., 2015. Therefore, stress tolerance QTLs for low and high temperature stress need to be stacked in the high yielding varieties for cultivation in the target environments. The inheritance of cold tolerance in rice is complex in nature and mechanism of the tolerance is difficult to explain using single or few genes. Gene mapping studies using various populations have detected role of many single genes and QTLs controlling cold stress tolerance during seedling stage of rice. More than 30 QTLs are reported to be involved in tolerance response for the trait in rice and they are present on different chromosomal locations (Kwak et al., 1984;Nagamine, 1991;Kim et al., 2000Kim et al., , 2014Misawa et al., 2000;Qian et al., 2000;Andaya and Mackill, 2003;Qu et al., 2003;Fujino et al., 2004;Zhan et al., 2005;Zhang et al., 2005;Andaya and Tai, 2006;Jiang et al., 2006Jiang et al., , 2008Han et al., 2007;Lou et al., 2007;Koseki et al., 2010;Wang et al., 2011;Suh et al., 2012). Pyramiding of these QTLs for cold tolerance breeding program is difficult due to more in numbers governing the tolerance response, lack of information on robust molecular markers and highly tolerant donors containing many QTLs for marker-assisted breeding. Screening and identification of highly tolerant donors for seedling stage cold tolerance from the germplasm pool is an important step for selection of donor parents in cold tolerance breeding program. Molecular markers, usually SSR or Simple Sequence Length Polymorphism (SSLP) are commonly being used for assessing rice genetic diversity (Panaud et al., 1996;Wu and Tanksley, 1996;Xiao et al., 1996;Olufowote et al., 1997;Thanh et al., 1999;Herrera et al., 2008;Pervaiz et al., 2010;Das et al., 2013;Babu et al., 2014;Anandan et al., 2016;Pradhan et al., 2016). Hence, genetic diversity available in parental lines should be assessed by using the trait linked SSRs and gene specific markers for seedling stage chilling stress tolerance before selection of parental lines. The genetic diversity and structure of the population for the trait of concern need to be studied for association mapping which could be useful in molecular breeding programs. For detecting a perfect marker-phenotype association, the population should not show spurious association or unequal relatedness within the population (Jagadish et al., 2007). Therefore, population structure (Q) with relative kinship (K) matrics analyses is used to check and correct the panel population composition for linkage disequilibrium (LD) mapping analyses (Yu et al., 2006). Thus, marker-based mixed linear model for better kinship estimate is considered appropriate for association mapping approaches that have shown to perform better than other model analysis.
Most of previous papers reported QTLs for the trait were derived through bi-parental mapping populations and mainly in japonica rice background. Bi-parental mapping is much resource consuming and very less number of alleles can be identified taking long time period with less resolution (Cardon and Bell, 2001;Flint-Garcia et al., 2003;Stich et al., 2005;Roy et al., 2006;Pradhan et al., 2016). These problems are eliminated or reduced by association mapping or linkage disequilibrium (LD) mapping approach. The purpose of LD mapping is to estimate the correlations between genotypes and desired phenotypic trait in a panel population containing genotypes from all the phenotypic classes. Association mapping, also known as linkage disequilibrium mapping, utilizes allelic variation in natural populations, and is capable of identifying many loci simultaneously for multiple traits (Flint-Garcia et al., 2005). The association mapping in rice has been reported for various traits like grain yield components (Agrama et al., 2007), agronomic traits (Huang et al., 2010), agro-morphologic and yield traits , tillering stage cold tolerance (Zhang et al., 2012), flowering and grain yield traits (Huang et al., 2012), grain quality traits (Zhao et al., 2013), drought tolerance traits (Muthukumar et al., 2015), salinity tolerance (Kumar et al., 2015), cold tolerance at germination and booting stages (Pan et al., 2015) mineral element in grain (Huang et al., 2015), early seedling vigor (Anandan et al., 2016), and high temperature stress tolerance . However, there is no such report till now for QTLs or genes for seedling stage cold tolerance through association analysis in indica rice. Therefore, in this study, a set of 304 shortlisted genotypes comprising mainly rice cultivars of high-elevation regions of India along with tolerant genotypes of other places of the country were analyzed for genetic diversity, population structure and association mapping using 60 linked markers for the trait. The robust markers showing strong effects and the identified highly tolerant donors could be the potential donor parents and markers in developing cultivars with seedling stage cold tolerance by molecular marker-assisted selection.

Plant Materials and Field Screening
The shorted materials comprising of 304 germplasm lines were pooled for the study from ICAR-National Rice Research Institute, Cuttack (NRRI) gene bank. Majority of the collection were from cold and higher elevation regions of the country including cold tolerant lines of NRRI (Table S1). The accessions were sown during December, 2014 in plots arranged in an augmented block design placing four check varieties (KalingaIII, Sahabhagidhan, Annada, and Vandana). A population size of three rows per germplasm line with a row length of 4 m and spacing of 20 × 15 cm were followed during sowing of seeds. Daily average minimum temperature was recorded during seedling stage of the crop. Scoring of germplasm  lines were done during the very low temperature period and recorded the scores for the change in leaf color, leaf rolling and plant growth. Screening of the germplasms for cold tolerance was done following IRRI-SES score (IRRI, 2014).

Screening under Control Condition
A panel of 66 genotypes representing all phenotypic classes were direct seeded in partially soil filled pots and uniform growth was maintained under green house chamber of RGAcum-Phytotron of the Institute for phenotyping chilling stress tolerance under controlled condition during June-July, 2015 ( Table 1). Genotypes were screened following the published protocol for chilling stress tolerance (Pradhan et al., 2015). As per the protocol, seedlings were grown in RGA-cum-Phytotron till three-leaf stage set at 25 • C temperature and photoperiod of 12 h. At this stage, the weak seedlings were removed and healthy seedlings were exposed to low temperature regime (LTR) treatments in the growth chamber with a setting of 75-85% relative humidity and 800 µ moles s −1 m −2 light intensity. The low temperature regime1 (LTR1) was by exposing the seedling to a temperature of 25 • C for 3 days; low temperature regime 2 (LTR2) with a gradual decrease in temperature to 15 • C from 25 • C and maintained at this temperature for 7 days; low temperature regime 3 (LTR3) with a gradual decrease in temperature to 8 • C from 15 • C and maintained at this temperature for 7 days; low temperature regime 4 (LTR4) with a gradual decrease in temperature to 4 • C from 8 • C and maintained at this temperature for 7 days; low temperature regime 5 (LTR5) with an exposure temperature of 4 • C for 14 days and low temperature regime 6 (LTR6) with an exposure temperature of 4 • C for 21 days. Complete randomized designs with two replications were followed for the experiment. The genotypes were evaluated using a modified IRRI-SES score under control facility. The modified score was given for field screening with miner modification as follows. 1 = Seedlings dark green; 3 = Seedling light green; 5 = Seedling with rolled leaves and light green to brownish yellow in color; 7 = Seedling with rolled leaves and brownish green in color; 9 = Seedling with rolled leaves with brown in color.

DNA Isolation and Selection of SSR Markers
Leaf samples were collected from 2 week old seedlings of 66 genotypes present in the panel population. Total genomic  DNA of all the studied genotypes were isolated using liquid nitrogen for grinding using CTAB extraction buffer (100 mM Tris-HCl pH 8, 20 mM EDTA pH 8, 1.3M NaCl, 2% CTAB) and chloroform-Isoamyl alcohol extraction followed by RNAase treatment and ethanol precipitation (Murray and Thompson, 1980). DNA concentration was estimated by using agarose gel electrophoresis. The sample DNA was diluted to ∼30 ng/µL. Earlier published reports on bi-parental mappings were used to select the seedling stage cold tolerance markers (58 linked SSR and 2 direct) and the distribution of markers covered all the chromosomes to illustrate the diversity.

PCR Amplification and Visualization of Markers Linked to Seedling Stage Chilling Stress
Polymerase chain reaction was performed by taking 20 µl aliquot using 1.5 mM Tris HCL (pH 8.75), 50 mM KCL, 2 mM MgCl 2 , 0.1% TrotonX-100, 200 µM each of dATP, dCTP, dTTP, dGTP, 4 pmole of each forward and reverse primers ( Table 2), 1 unit of Taq polymerase and 30 ng of genomic DNA. A Programmable Thermal Cycler was used for amplification of genomic DNA samples (Veriti, Applied BioSciences). First, the reaction mixture was denatured for 4 min at 94 • C and then continued to 35 cycles of 1 min denaturation at 94 • C, 1 min annealing at 55 • C, 1 min extension at 72 • C, and then a final extension for 10 min at 72 • C. Agarose gel of 2.5% containing 0.8 µg/ml Ethidium Bromide was used for electrophoresis. Aliquots of 10 µl of the products from PCR amplification were loaded in the agarose gel. Size of amplicons was determined by using 50 bp DNA ladder. The gel was run at 60 volts (2.5 V/cm) in 1X TBE (pH 8.0) for 4 h and photographed using a Gel-Doc System (SynGene).

Genetic Diversity, Population Structure, and Linkage Disequilibrium Mapping
Data were scored on the basis of presence or absence of the alleles for each genotype-primer combination and arranged in a binary data matrix as discrete variables. PowerMarker Ver3.25 program was used for data analysis to generate number of alleles, allele frequency, gene diversity, heterozygosis, and polymorphic information index (PIC; Lu et al., 2005). STRUCTURE 2.3.4 software a model based approach was used for data analysis to obtain possible population structure (Pritchard et al., 2000). We ran STRUCTURE 2.3.4 software with model parameter set of "possibility of admixture and allele frequency correlated" with a burn-in period of the 150,000 followed by 150,000 Markov Chain Monte Carlo (MCMC) replications. Each K-value was run for 10 times with K-value ranging from 1 to 10. The optimum K-value was determined by plotting the log posterior probability data to the given K-value. The maximal value of L (K) was identified using the exact number of sub-populations. The model choice criterion to detect the most probable value of K was K, an ad-hoc quantity related to the second-order change of the log probability of data with respect to the number of clusters inferred by STRUCTURE (Evanno et al., 2005). Structure Harvester was used for estimation of the K-value as function of K showing a clear peak as the optimal K-value (Earl and Von, 2012). A total of 60 linked markers distributed on 12 chromosomes were adopted to divide the accessions into different groups with the membership probabilities threshold of 0.80 as well as the maximum membership probability among groups. Those accessions with <0.80 membership probabilities were retained in the admixed group. NEI coefficient for dissimilarity index was calculated by constructing unweighted neighbor joining un-rooted tree (Nei, 1972) with bootstrap value of 1,000 by using DARwin5 software (Perrier and Jacquemoud-Collet, 2006). The Principal Coordinate analysis of the germplasm lines was performed as per the standard published method followed in published paper . GenAlEx 6.5 software was used to assess the presence of molecular variance within and between the population structures using analysis of molecular variance (AMOVA; Peakall and Smouse, 2012). F statistics including deviations from Hardy-Weinberg expectation across the whole population (F IT ), deviation from Hardy-Weinberg expectation within a population (F IS ) and correlation of alleles between subpopulation (F ST ) was calculated. The hypothesis of the association of SSR markers with seedling stage cold tolerance was tested using a general linear model (GLM) and mixed linear model (MLM) in the program TASSEL 5 (Bradbury et al., 2007). Linkage disequilibrium plot was constructed using LD measured r2, between pair of markers is plotted against the FIGURE 2 | Genotype-by-trait biplot analysis of 66 genotypes for first two principal components. LTR1-low temperature regime of 25 • C for 3 days; LTR2-low temperature regime of 15 • C for 7 days; LTR3-low temperature regime of 8 • C for 7 days; LTR4-low temperature regime of 4 • C for 7 days; LTR5-low temperature regime of 4 • C for 14 days; LTR6-low temperature regime of 4 • C for 21 days. The numbers in the figure represent the serial number of the genotypes enlisted in Table 1.
distance between the pair. FDR adjusted p-values (q-values) were estimated using statistical software SPSS version 16.

Phenotypic Screening of Rice Genotypes for Seedling Stage Cold Tolerance under Field Situation
Phenotypic screening experiment was performed using 304 shortlisted rice germplasm lines during 2014-15 dry season at ICAR-National Rice Research Institute (NRRI), Cuttack, Odisha, India. The shortlisted genotypes were pooled on the basis of its area of adaption in high-altitude area of the country and cold tolerance reaction of the lines from NRRI germplasm catalog. These genotypes were raised during dry season to coincide with the seedling stage of the crop to the coldest period of the year at the Institute site (Figure 1). The field screening results exhibited 86 genotypes to be tolerant to seedling stage cold stress with a score of 1-5. The screening score ranged from SES score of 1 (Geetanjali, KalingaIII, Langma and Umleng-1) to score 9 (Aujari, CR Dhan 505, CR 2340-1, CR 3820-4-5-3-4-1, 9022, Sarpung, Kala Joha, Champa, Hema, Sadabahar, Swarnaprava, Mayang khang, Mayang khang-2, and Yentie) (Table S1). On the basis of subtle changes in leaf color, the shortlisted genotypes for seedling stage cold tolerance under field condition showed 1.3% to be highly tolerant (Score 0-1); 11% tolerant (score 2-3); 16% moderately tolerant (score 4-5); 67% sensitive (score 6-7), and 5% very highly sensitive (score 8-9). A panel population containing 66 genotypes was constituted representing all phenotypic groups for phenotyping under control facility and association mapping of the trait.  Table 1.

Phenotyping of Panel Population for Seedling Stage Cold Tolerance under Control Facility
Under field screening, 86 genotypes were found to be tolerant to seedling stage cold with a score of 1-5. Rest, 218 genotypes were observed to be susceptible to cold stress. The screening experiment resulted from control facility by using a cold chamber is presented in Table 1. Under LTR1 exposure, susceptible check variety Sahabhagidhan became brownish yellow with rolled leaf showing SES score of >5 whereas rest 65 genotypes showed normal growth and color. The tolerant genotypes under LTR1 were further exposed to LTR2 to separate out susceptible genotypes in this temperature regime. The 15 separated genotypes were Tejaswini Changlei-1, Manipur local 2, Phouurel, Photem, Manipur local1, Napdai, Sadabahar, Swarnaprava, Swarna, Pratikhya, Tapaswini, Gayatri, Ranidhan, and Naveen showed rolled leaf, yellow to brown in color and reduced growth. Hence, the above two types of susceptible genotypes were categorized into two susceptible sub-groups.
Genotypes showing >5 SES score at evaluation temperature of >15 • C for 7 days as highly susceptible (HS) types while genotypes exhibiting >5 SES score at above chilling stress (>8 • C) and <15 • C for 7 days as moderately susceptible (MS) genotypes. The remaining 49 genotypes showing a score of 1-5 with regard to leaf color and leaf rolling at 8 • C were further exposed to LTR3, subsequently under which 20 genotypes were observed with score of >5. These 20 genotypes can be grouped as moderately tolerant group to chilling stress at seedling stage ( Table 1). Rest 29 genotypes with a score of 1-5 were further exposed to LTR4 (chilling stress up to 14 days at 4 • C), only 2 genotypes namely Tabadugu and Krishnahamsa were found to be with SES score of >5 for chilling stress tolerance under this temperature regime. Hence, Tabadugu and Krishnahamsa were categorized as a different group in relation to previous classes and termed as tolerant group. Rest 27 genotypes with 1-5 score were further exposed to 21 days under chilling stress for further grouping under tolerant sub-group (LTR5). At 4 • C for 21 days, genotypes exhibiting >5 SES score were Kalinga-3, Radang, Govind, Rungpchi, Ajay, Kamesh, Paun, Manipuridhan, MR37, Umleng-2, Jamak, Charmui, Umbo, Langme-1, Mopu, Tazek, Serum, Bamak, Lagmin, Itanagardhan, and Langme-2. This sub group can be categorized into highly tolerant (HT) group. Genotypes Geetanjali, Langma, and Umleng-1 exhibited a score   Figure S1.

Genetic Relatedness by Biplot, Principal Coordinate, and Cluster Analyses
The phenotyping data of the panel population for chilling stress tolerance were used to produce genotype-by-trait biplot graph for depicting the genotypes in the first two principal components (Figure 2). The first principal component showed 82.51% of variation, while second component explained 11.03% of the total variability. Among the temperature regime studied, LTR3 contributed maximum toward diversity, followed by LTR5 and LTR4 (Figure 2). The top left (Ist quadrant) and bottom left (4th quadrant) accommodated 24 genotypes, which were tolerant to highly tolerant in response to chilling stress during seedling stage. The two tolerant genotypes in 4th quadrant were clearly separated out from the seven highly tolerant genotypes. The 15 genotypes in the 1st quadrant showed better tolerance than the seven highly tolerant genotypes in 4th quadrant. The 2nd quarter contained all the intolerant genotypes to the stress exhibiting higher bronzing symptoms along with few moderately tolerant genotypes. These moderately tolerant genotypes were placed near the axis away from the intolerant ones. The 3rd quadrant genotypes that constituted only moderately tolerant genotypes were marginally better tolerant than the genotypes in 2nd quadrant. The encircled area in the Figure 2 exhibited 22 highly tolerant genotypes. Principal coordinate analysis (PCoA) was performed using the linked markers for determining the genetic relatedness among the genotypes (Figure 3). In PCoA, the desirable genotypes are depicted in the encircled area most of which are located in 1st quadrant. The pattern of distribution of genotypes showed a clear grouping based on tolerance level to chilling stress. The genotypes Umleng-2, Mopu, Umleng-1, Jamak, Manipuri dhan, Tazek, Langma, Paun, MR37 (Farmer selection), Rungpchi, Lagmin, Langme-2, Itanagar dhan, and Tazek clustered together depicted in the encircled area were high to VHT to chilling stress (Figure 3). The first two axis of differentiation explained 19.85 and 14.11% of the total variation, respectively.
Cluster analysis was carried out to assess genetic distance and the dissimilarity matrix-using UPGMA method. The tree constituted of three major clusters (Figure 4). While the phenotyping results of panel population was compared with the grouping pattern, significant correlation was observed ( Figure 4A). Cluster 1 possessed 19 genotypes accommodating high to VHT genotypes (pink color) in the tree except Tabadugu (Figure 4A). Similarly, the cluster 2 contained all the sensitive genotypes to chilling stress except IR-64 and Virendra. The moderately tolerant and tolerant germplasm lines were mainly observed in cluster 3 along with few highly tolerant and susceptible types indicating grouping of the admixture type landraces for the trait. However, the UPGMA tree correlated with population structure analysis, the entire genotypes were clearly categorized into three groups and depicted in three colors ( Figure 4B). The first cluster consisted all the 19 genotypes belonging to sub-population 1 (red colors). Similarly, cluster 2 (blue color) consisted 12 genotyped classified as SP2 by structure analysis, whereas cluster 3 (green) consisted rest 35 genotypes that were classified as SP3.

Genetic Diversity
The panel containing 66 genotypes, comprising tolerant and susceptible types were genotyped using 58 SSR linked and 2 gene specific markers for the trait. All the loci used for genotyping the panel rice germplasm lines and their genetic diversity parameters obtained are presented in Table 3. The total numbers of amplicons were 222 by using 58 co-dominant and 2 dominant markers. An average of 3.6667 alleles per locus was detected with a range of 1-10 per marker with the highest number of 10 alleles from RM 1812 in the panel for seedling stage chilling stress tolerance. The mean PIC value was observed to be 0.4540 with minimum value of 0.000 (RM200, RM315, RM5221, RM6091, RM6947, and RM 6547) and maximum of 0.7787 (RM1812). The observed average heterozygosity (H o ) was 0.1912 which ranged between 0.00 and 1.0. It is observed that 42 markers exhibited the level of H o more than zero, while remaining 18 showed zero values. The average heterozygosis or average gene diversity (H e ) was found to be 0.5069 which varied from 0.0000 (RM319, RM522, RM315, RM200, RM6091, RM6947, and RM6547) to 0.7990 (RM1812). The major allele frequency of these chilling stress tolerance linked polymorphic markers ranged from 0.2879 to 1.0000 with an average of 0.5792 (Table 3).

Population Structure
The panel population was divided into three sub-populations for chilling stress tolerance on the basis of analysis by STRUCTURE software (Figures 5B, 6). The population panel was analyzed for genetic structure on the Bayesian clustering approach taking probable sub-populations (K) and selecting higher K-value, an ad-hoc quantity related to the second order change of the log probability of data for the number of clusters detected by Structure (Evanno et al., 2005). A high K peak value of 95.6 was observed among the assumed K at K = 3 as per the Evano table output (Figure 5A). The sub-population 1 (SP1) contained 18 genotypes with 14 pure and 4 admixture types, accommodating highly and VHT genotypes to seedling stage cold tolerance. A total of 12 genotypes present in sub-population 2 (SP2) representing susceptible sub-population were with all pure type to the sub-population (Figures 5B, 6 and Table 4). The third sub-population can be grouped as moderately tolerant group with 24 tolerant, 6 high to VHT, and 6 susceptible germplasm lines in it. This sub-population has 19 pure and 17 admixture ones genotypes in it. Maximum allele frequency divergence between populations was observed in SP1 and SP2 (0.1664) while within divergence sub population was highest in SP3 (0.242). The fixation index values (F ST ) of the sub-populations were found to be 0.4947, 0.4663, and 0.3297 for SP1, SP2, and SP3, respectively.    Figure S1).

Analysis of Molecular Variance (AMOVA)
The three populations obtained through structure analysis were used for genetic variation between and within the clusters using AMOVA (  Figure  S1).

Association of Marker Alleles with Cold Tolerance
The marker-trait associations for cold tolerance with various temperature regime and treatment duration were calculated using GLM and MLM (Q+K) model of TASSLE5 software. The comparisons were filtered with p < 0.05. The r 2 -values varied from 0.0594 to 0.4062 with an average of 0.164 by using GLM, whereas the average reduced to 0.0919 with upper border value of 0.283 and lower border 0.062 by using MLM analysis ( Table 7 and Table S2). Among 60 markers used, 48 markers were associated with different level of cold tolerance by using GLM, whereas the number reduced to 23 with MLM model at p < 0.05 and r 2 > 0.05 (Table 6). In total, 130 comparisons with GLM and 8 comparisons with MLM were significant with r 2 > 0.10 at p < 0.05. All 130 comparisons significant with GLM at p < 0.05 FIGURE 5 | (A) Graph of estimated membership probability fraction for K = 3. The maximum of ad-hoc measure K determined by structure harvester was found to be K = 3, which suggested that the entire population can be categorized into three sub-groups (SP1, SP2, and SP3) and (B) population structure of a panel possessing 66 genotypes based inferred ancestry using 60 molecular markers (K = 3).
were also significant at p < 0.01, but only 5 comparisons were significant with MLM at p < 0.05 and r 2 > 0.10. Cold tolerance at different temperatures like 25, 15, 8, and 4 • C at 7, 14, and 21 days treatment were associated with the marker data. Eight markers namely, RM152, RM341, RM50, RM4154, RM245, RM13335, RM282, and RM1341 were associated with cold tolerance at 4 • C for 21 days at p < 0.01 and r 2 > 0.10 with GLM analysis. Similarly, 28, 30, 34, and 30 numbers of markers were associated with tolerance to 15 • C for 7 days, 8 • C for 7 days, 4 • C for 7 days, and 4 • C for 14 days, respectively ( Table 6). Higher F-value and lower p-value with high r 2 indicated the positive association with the trait. Further, MLM analysis was performed to achieve more precise association, considering the kinship value. This showed a strong marker-trait association with phenotypic variance of 11.31% among RM558 and tolerance for 7 days at 15 • C to 28.33% among RM1341 and tolerance for 21 days at 4 • C considering p < 0.01.
TASSEL analysis also evidenced association of some markers with all the studied temperature regimes and duration, whereas some other markers were either regime or duration specific. The markers like RM1347, RM328, RM152, RM341, RM50, RM2634, RM4112, RM5310, RM7179, RM3701, RM104, RM9, and RM1211 were positively associated with all treatments considered, except at 25 • C ( Table 7 and Table S2). Twelve markers namely, RM3375A, RM5746, RM286, RM84, RM561, RM253, RM284, RM239, RM256, RM1812, RM558, and RM173 were associated with tolerance at 15 • C for 7 days to 4 • C for 14 days. The markers RM245, RM3602, RM493, RM13335, RM282, and RM5704 did not show any association with tolerance for 7 days at 15 • C but all other treatments. Three temperature regime specific markers RM297, RM13341 and RM506 had been detected which were associated with all treatment durations at 4 • C, whereas IN11D1 and RM590 were specific for 15 • C only. Further, RM6651 and RM22491 were associated with 15 • C and 8 • C only. The markers IN1C3, RM1113, RM3648 and RM85 were associated with 8 • C for 7 days and 4 • C for 7 and 14 days but not 21 days. Two markers RM7003 and RM305 showed association with only 21 days treatment at 4 • C. However, no significant association was observed for 3 days treatment at 25 • C. The QQ plot also confirmed significant association of markers for all temperature regime and treatment duration except 25 • C treatment (Figure 7). The linkage disequilibrium decay plot for seedling stage cold tolerance has been plotted using r 2 -value between pair of markers and distance between the pair ( Figure S2).

DISCUSSION
Low temperature stress during seedling stage is one of the serious yield reducing abiotic factors in dry season rice, particularly boro rice of India and cultivation in the high altitude areas like hill rice. The effect of seedling stage cold prolongs the growth duration of the rice plant, subsequently delay in flowering period of boro and dry season rice that coincides with high temperature stress period resulting in drastic reduction in yield . Therefore, cold tolerance breeding for seedling stage is important in boro and dry season rice. Results of field screening experiment followed by controlled screening method under RGA-cum-Phytotron, we could identify 50 germplasm lines with seedling stage cold tolerance. The genotypes were classified into six classes based on their tolerance to cold stress after exposure to six temperature regimes. Clustering by using results on molecular markers and genotype-trait biplot analysis exhibited grouping of genotypes basing on their tolerance to seedling stage cold tolerance (Figures 2-4). The tolerant genotypes were grouped into many sub-groups depending upon their level of tolerance, which might be the result of expression of different gene(s)/QTL(s) for seedling stage cold tolerance and their possible presence in the studied materials. Hence, the population structure of the panel population for the trait is most important. Earlier studies also indicated screening and identification of cold tolerant genotypes in rice (Kwak et al., 1984;Nagamine, 1991;Kim et al., 2000Kim et al., , 2014Misawa et al., 2000;Qian et al., 2000;Andaya and Mackill, 2003;Qu et al., 2003;Fujino et al., 2004;Zhan et al., 2005;Zhang et al., 2005;Andaya and Tai, 2006;Jiang et al., 2006Jiang et al., , 2008Han et al., 2007;Lou et al., 2007;Koseki et al., 2010;Wang et al., 2011;Suh et al., 2012;Pradhan et al., 2015). Germplasm lines like AC 43281 (Langma), AC 43291 (Umleng 1), and Geetanjali were observed to be VHT to seedling stage cold tolerance. These results are also confirming the earlier screening results for chilling stress tolerance in rice (Pradhan et al., 2015).
The principal component and coordinate analysis placed the tolerant and non-tolerant genotypes into various spots in the four quadrants. The location of the genotypes with respect to origin provides clue for tolerance to seedling stage cold tolerance. Thus, distribution of genotypes in different quadrant confirmed the presence of variation for the trait. The UPGMA tree, categorized the highly tolerant and tolerant classes on the basis of banding pattern of 60 cold stress linked molecular markers. The genotypes could be categorized into separate clusters which are in line with the seedling stage cold tolerance phenotype groups. Various groups and sub-groups of the genotypes obtained based on banding pattern of cold tolerance linked markers suggested the presence of many genes/QTLs in the panel population. In our study, a moderate level of genetic diversity was detected for seedling stage cold tolerance. Our results on the level of genetic diversity for the trait is similar to the other results of moderate diversity parameters reported earlier for various traits (Agrama and Eizenga, 2008;Jin et al., 2010;Chen et al., 2011;Zhang et al., 2011;Shah et al., 2013;Singh et al., 2013). However, very rich genetic diversity values for agro-morphologic traits were also reported earlier in rice (Garris et al., 2003;Zhao et al., 2013;Salgotra et al., 2015).
We demonstrated the appropriateness of the suggested panel population for association mapping and kinship study basing on population structure and relatedness for seedling stage cold tolerance. The phenotyping and genotyping of the panel population using 60 linked markers for the trait clearly categorized the study materials into different groups, suggesting their differential response to cold stress tolerance (Figures 2-4). This heterogeneity favored the presence of linkage disequilibrium and increased the chance of recovering marker-trait association. Marker-phenotypic trait association was also detected earlier  showing the potential value of germplasm in heterogeneous collections (Gebhardt et al., 2004;Lu et al., 2005;Caicedo et al., 2007;Zhang et al., 2009Zhang et al., , 2011Zhao et al., 2013;Anandan et al., 2016;Pradhan et al., 2016). Boro and dry season rice affected by low temperature stress during seedling need attention for incorporating seedling stage cold tolerance. The application of association mapping results for this trait will be helpful in selecting the molecular markers to be used in marker-assisted breeding of the multiple QTLs/genes responsible for the trait. Therefore, different landraces with differential response to the stress are needed for association mapping of the trait. In our results, we detected a lower value of alpha (α = 0.1079) from which we can infer that in most of the landraces, the trait had a common primary ancestor with few admix individuals in each sub-population. The inferred ancestry indicated that small effects QTLs for the trait present in different landraces might be pooled together from long time ago through intercrossing of many landraces naturally as a result of which few landraces possess many QTLs that are strongly tolerant to the stress. Similar views have also been provided by previous reports (Mather et al., 2004;Zhao et al., 2013;Pradhan et al., 2016). Values of F ST were very high in SP1 and SP2 (0.403), SP3 and SP1 (0.281) and SP3 and SP2 (0.264) when combination of inter sub-populations are taken together, thus suggesting higher genetic differences between germplasm accessions. It has been demonstrated that populations and individuals with higher F ST values produce better variable materials when combined with lines from different genetic diversity estimates (Watkins et al., 2003). The significant F ST among the clusters indicate a real variation in these clusters, and attempt to pyramid the QTLs governing the trait may improve further tolerance for seedling stage cold tolerance. Similar suggestion was also provided by earlier workers for increasing heterosis for grain yield in rice (N'Goran et al., 2000). Cold stress exposure at different temperature regimes of the tested panel genotypes was observed to be associated with the marker data. Higher F and lower p-value with high r 2 were detected for 28, 30, 34, 30, and 8 markers associated with tolerance to LTR 2 (15 • C for 7 days), LTR 3 (8 • C for 7 days), LTR 4 (4 • C for 7 days), LTR 5 (4 • C for 14 days), and LTR 6 (4 • C for 21 days), respectively. This indicated more numbers of markers detected by TASSEL using both GLM and MLM analyses providing a robust marker-phenotype association in the present study which was also evident from Q-Q plot (Figure 7). This also suggests that, the observed marker-trait association resulted possibly from multiple introgressions or accumulation of QTLs from the landraces by natural hybridization over a long time period ultimately exhibiting a higher level of tolerance. Similar results on multiple introgressions of tolerant QTLs for high temperature stress from the landraces were also earlier described in rice . In this study, we found a strong marker-trait association by both MLM and GLM models of TASSEL analysis in all the studied temperature regimes with some regime or duration specific markers association except LTR1. This confirms the effectiveness of the linked markers for cold tolerance not showing any significant association with temperature treatment at 25 • C. Twelve markers namely, RM3375A, RM5746, RM286, RM84, RM561, RM253, RM284, RM239, RM256, RM1812, RM558, and RM173 were associated with tolerance at 15 • C for 7 days to 4 • C for 14 days. Tolerance in these temperature regimes suggests that the markers associated at these low temperature exposures were with QTLs/genes conferring tolerance from mild to highly tolerance response to the stress. The phenotype groups associated for these primers were MS, MT, and T types of phenotypes. Similarly, markers like RM1347, RM328, RM152, RM341, RM50, RM2634, RM4112, RM5310, RM7179, RM3701, RM104, RM9, and RM1211 were positively associated with all treatments considered except 25 • C. This suggests that tolerance to all the temperature regimes means germplasm lines possessing many tolerance conferring QTLs for cold tolerance. Here, these primers are associated with four phenotype groups observed for the trait. This is also evidenced from the earlier bi-parental mapping population indicating role of QTLs like qCTS11.1, qCTS9, qCTS6-1, qSCT2, FIGURE 7 | Quantile-Quantile (QQ) plot and distribution of marker-trait association from MLM analysis at different low temperature regimes.
qSCT1a, qCTS-3.1, qCTS-2, qCTS12.1, and qCTS-1b for seedling stage cold tolerance (Andaya and Mackill, 2003;Long-zhi et al., 2004;Lou et al., 2007;Wang et al., 2011;Kim et al., 2014). Similarly, markers RM245, RM3602, RM493, RM13335, RM282, and RM5704 did not show any association with tolerance for 7 days at 15 • C but all other treatments. This indicates the association of tolerant phenotypes excluding the chilling sensitive alleles. These associations of primers with tolerance phenotypes indicate that the tolerant genotypes chosen here were from all tolerant groups and may possess many tolerant QTLs for the trait. Similar results on high temperature stress tolerance governed by multiple tolerant QTLs were reported earlier in rice . Previous research results on cold stress which suggested presence of many cold tolerant QTLs on chromosome 1, 2, 3, 4, 5, 6, 8, 9, 10, 11, and 12 obtained by using various mapping populations (Andaya and Mackill, 2003;Long-zhi et al., 2004;Lou et al., 2007;Wang et al., 2011;Suh et al., 2012;Kim et al., 2014;Bonnecarrere et al., 2015;Pradhan et al., 2015). Also, our PCA, PCoA, and UPGMA trees revealed presence of highly tolerant lines in the encircled area (Figures 2-4). The study showed that out of the earlier reported QTLs (Table 1), only nine QTLs namely, qCTS11.1, qCTS9, qCTS6-1, qSCT2, qSCT1a, qCTS-3.1, qCTS-2, qCTS12.1, and qCTS-1b were significantly involved for broad level of cold tolerance at seedling stage, whereas the other QTLs were less effective. This may be due to the linked markers chosen for our study which were earlier reported for various related traits through bi-parental mapping approach. When the association mapping was studied in multiple genotypes, it revealed robustness of these QTLs and linked markers suggesting their further use in marker-assisted seedling stage cold tolerance breeding program.

CONCLUSION
From the experiment, a moderate level of genetic diversity for seedling stage chilling tolerance was noticed by using the panel population. Results of the STRUCTURE analysis revealed that the entire population could be grouped into three subpopulations and further detected that most of the landraces had a common primary ancestor with few admix individuals for the trait. The donor lines in the panel exhibited the presence of different QTLs representing like whole genome diversity for the expression of tolerance. The significantly associated markers like RM1347, RM328, RM152, RM341, RM50, RM2634, RM4112, RM5310, RM7179, RM3701, RM104, RM9, RM1211, RM245, RM3602, RM493, RM1335, RM282, and RM5704 were significantly associated at chilling stress of 8 • C to 4 • C for 7-21 days duration. Thus, the primers linked to the seedling stage cold tolerance QTLs namely qCTS9, qCTS-2, qCTS6.1, qSCT2, qSCT11, qSCT1a, qCTS-3.1, qCTS11.1, qCTS12.1, qCTS-1b, and CTB2 need to be pyramided for development of strongly chilling tolerant variety.