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).

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 RGA-cum-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⁻¹ m⁻² 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.

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.

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₂, 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).

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 (Pradhan et al., 2016). 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 distance between the pair. FDR adjusted p-values (q-values) were estimated using statistical software SPSS version 16.

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.

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 of 1–5 after 21 days at chilling stress which are categorized as very highly tolerant (VHT) group. The highly tolerant genotypes like Kalinga III (positive check) could survive for 20 days and VHT genotypes AC 43281 (Langma), AC 43291 (Umleng 1) and Geetanjali even survived for 24 days with a score of 7 under exposure to LTR6 chilling stress. Representative figures of different classes of tolerant and susceptible genotypes are depicted in the Figure S1.

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.

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).

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. Further, the program exhibited a lower value of alpha (α = 0.1079) for the population panel. The distribution pattern of α-value in the panel population and distribution of F_ST values in the sub-populations are shown in Supplementary File (Figure S1).

The three populations obtained through structure analysis were used for genetic variation between and within the clusters using AMOVA (Table 5). The analysis accounted 30% of the variation among populations, 70% among individuals while no variation detected within individuals in the panel population. Hardy-Weinberg deviation in the population was detected using Wright's F statistic. Similar value of F_IS and F_IT for all the 60 markers loci were 1.00, while F_ST was 0.298 among populations. Pair wise F_ST values showed significant differentiation among all the pairs of sub-populations ranging from 0.001 to 0.403 indicating that all the three groups were significantly different from each other. A higher F_ST values were observed for pair wise sub-populations in SP1 and SP2 (0.403), SP3 and SP1 (0.281) and Sp3 and SP2 (0.264) when three inter sub-populations are considered. The F_ST values and their distribution pattern show clear differentiation of sub populations from each other (Figure S1).

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²-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² > 0.05 (Table 6). In total, 130 comparisons with GLM and 8 comparisons with MLM were significant with r² > 0.10 at p < 0.05. All 130 comparisons significant with GLM at p < 0.05 were also significant at p < 0.01, but only 5 comparisons were significant with MLM at p < 0.05 and r² > 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² > 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² 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²-value between pair of markers and distance between the pair (Figure S2).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.