Phenotypic and Molecular Characterization of Rice Genotypes’ Tolerance to Cold Stress at the Seedling Stage

: Rice plants are affected by low-temperature stress during germination, vegetative growth, and reproductive stages. Thirty-nine rice genotypes including 36 near-isogenic lines (NILs) of BRRI dhan29 were evaluated to investigate the level of cold tolerance under artiﬁcially induced low temperature at the seedling stage. Three cold-related traits, leaf discolouration (LD), survivability, and recovery rate, were measured to determine the level of cold tolerance. Highly signiﬁcant variation among the genotypes was observed for LD, survivability, and recovery rate. Three NILs, IR90688-74-1-1-1-1-1, IR90688-81-1-1-1-1-1, and IR90688-103-1-1-1-1-1, showed tolerance in all three traits, while IR90688-118-1-1-1-1-1 showed cold tolerance with LD and recovery rate. IR90688-92-1-1-1-1-1, IR90688-125-1-1-1-1-1, IR90688-104-1-1-1-1-1, IR90688-124-1-1-1-1-P2, IR90688-15-1-1-1-1-1, and IR90688-27-1-1-1-1-1 showed signiﬁcantly higher yield coupled with short growth duration and good grain quality. Genetic analysis with SSRs markers revealed that the high-yielding NILs were genetically 67% similar to BRRI dhan28 and possessed cold tolerance at the seedling stage. These cold-tolerant NILs could be used as potential resources to broaden the genetic base of the breeding germplasm to develop high-yielding cold-tolerant rice varieties.


Introduction
Compared to other cereal crops, rice is more sensitive to low-temperature stress (LTS) as it has originated from tropical regions [1][2][3][4]. Cold stress negatively affects rice plants throughout various growth stages, from germination to maturity, and causes significant yield losses because of poor germination and seedling establishment, stunted growth pattern, non-vigorous plants, vast spikelet sterility, delay in flowering, and lower grain filling in the temperate, subtropical, and even in the tropical rice-growing regions [5][6][7][8][9][10]. Boro rice, which is grown from November to May in Bangladesh and some parts of Eastern India, is affected by LTS at both the seedling and reproductive stage. In the northern and north-eastern districts of Bangladesh, LTS at the seedling stage is a major problem for rice cultivation that seriously affects crop establishment. Seedling mortality occurs 10-90% of the time due to severe cold injury during transplanting in late December to

Sources of Materials Used
Thirty-six BC 3 F 5 near-isogenic lines (NILs) derived from a cross between a Korean japonica variety (Jinbubyeo) and a Bangladeshi indica-type high-yielding rice variety (BRRI dhan29), along with four check varieties (BR1, BR18, BRRI dhan28, and BRRI dhan29), were evaluated for seedling stage cold tolerance and agronomic performance. The same set of NILs, along with their parents and cold-tolerant and susceptible check varieties, were used in genotyping to explore the genetic variability/similarity that exists among them. The seeds of the NILs and the check varieties (BR1, BR18, BRRI dhan28, and BRRI dhan29) were collected from the International Rice Research Institute (IRRI) and the Bangladesh Rice Research Institute (BRRI), respectively.

Evaluation of Near Isogenic Lines (NILs) for Cold Tolerance at Seedling Stage
The cold screening was performed using a cold water tank in the cold screening laboratory of the Bangladesh Rice Research Institute (BRRI) following the protocol described in Khatun et al. [32] (Figure 1). Briefly, 10 pre-germinated seeds of each genotype were sown in one-row plots spaced at 3.0 cm in plastic trays of 60 cm × 30 cm × 2.5 cm in size. The trays were filled with gravel and crop residue-free granular and fertilized soil. After seeding, a thin layer of fine granular soil was used to cover the germinating seeds. The seedlings were allowed to grow at ambient temperature. At the three-leaf stage (approximately 12-15 days after seeding depending on the ambient temperature), the plastic trays were placed in a cold water tank pre-set at 13 °C. Sensitivity to cold stress was recorded using arbitrary leaf discolouration (LD) scores (1 to 9), survivability, and recovery rate as described in Biswas et al. [4]. The LD score and survivability were recorded at seven days of cold treatment or when the susceptible check variety (BR1) died. After LD and survivability recording, the trays were placed at ambient temperature in a sunny place. The recovery rate was recorded 7 days after the removal of cold stress. The survivability and recovery rate were calculated as the percentage of the green plants to the total number of plants used in the screening activity following the formula given below: Survivability rate (%) = No. of Green plants × 100 No. of plants treateḋ (1) Recovery rate (%) = No. of recovered plants × 100 No. of plants treated (2) The analysis of variance (ANOVA) and mean comparison of the LD scores, survivability, and recovery rate were performed using the analytical software Statistix version 10. Spearman correlation coefficients between the traits were calculated using the Corplot 0.92 package in R 3.6.2.

Evaluation of Near Isogenic Lines (NILs) for Agronomic Performance
All 36 NILs along with four check varieties were grown in the field following augmented RCB design [33]. At 35 days old, seedlings were transplanted with a single seedling per hill at a spacing of 20 cm × 20 cm in 5.4 m 2 plots. The standard crop management protocol described in BRRI [34] was followed throughout the crop growth period. Observations on different morphological characters viz. plant height (PH), number of effective tillers per hill (ETPH), panicle length, number of fertile grains per panicle (FGPP), spikelet After seeding, a thin layer of fine granular soil was used to cover the germinating seeds. The seedlings were allowed to grow at ambient temperature. At the three-leaf stage (approximately 12-15 days after seeding depending on the ambient temperature), the plastic trays were placed in a cold water tank pre-set at 13 • C. Sensitivity to cold stress was recorded using arbitrary leaf discolouration (LD) scores (1 to 9), survivability, and recovery rate as described in Biswas et al. [4]. The LD score and survivability were recorded at seven days of cold treatment or when the susceptible check variety (BR1) died. After LD and survivability recording, the trays were placed at ambient temperature in a sunny place. The recovery rate was recorded 7 days after the removal of cold stress. The survivability and recovery rate were calculated as the percentage of the green plants to the total number of plants used in the screening activity following the formula given below: . Survivability rate (%) = No. of Green plants × 100 No. of plants treated (1) Recovery rate (%) = No. of recovered plants × 100 No. of plants treated (2) The analysis of variance (ANOVA) and mean comparison of the LD scores, survivability, and recovery rate were performed using the analytical software Statistix version 10. Spearman correlation coefficients between the traits were calculated using the Corplot 0.92 package in R 3.6.2.

Evaluation of Near Isogenic Lines (NILs) for Agronomic Performance
All 36 NILs along with four check varieties were grown in the field following augmented RCB design [33]. At 35 days old, seedlings were transplanted with a single seedling per hill at a spacing of 20 cm × 20 cm in 5.4 m 2 plots. The standard crop management protocol described in BRRI [34] was followed throughout the crop growth period. Observations on different morphological characters viz. plant height (PH), number of effective tillers per hill (ETPH), panicle length, number of fertile grains per panicle (FGPP), spikelet fertility (%), thousand-grain weight (g), growth duration (days), and grain yield per hill (g) were recorded at the appropriate growth stage of the rice plant following the descriptor for rice [35]. The data collected for all eight quantitative characteristics were subjected to analysis of variance for augmented block design according to the method suggested by Federer and Raghavarao [36]. The major descriptive statistics such as mean, range, and coefficient of variation were estimated following Panse and Sukhatme [37]. The plot means were analysed using standard statistical analysis suggested by Federer [33] and elaborated by Sharma [38]. The means of the genotypes were adjusted with the block effects, which were measured from the replicated check plots following the formula: where, Rj = effect of jth block, Bj = mean of all check-in the jth block, and M = grand mean of the check.
Principal component analysis (PCA) was performed using the R package FactoMineR and visualized in biplot using the ggplot2 package.

Genetic Analysis of Near-Isogenic Lines (NILs) of BRRI dhan29 Using SSR Markers
A total of 58 simple sequence repeat (SSR) markers (Table S1) randomly distributed over 12 chromosomes of rice were used in the genotyping of 42 genotypes, including 36 NILs and their parents (BRRI dhan29 and Jinbubyeo), one susceptible check variety (BR1), and two tolerant check varieties (BR18 and Hbj. B. VI). Genomic DNA was isolated from young leaf tissues following the modified mini-scale method as described by Syed et al. [39]. The polymerase chain reaction (PCR) and gel electrophoresis were performed as described by Syed et al. [40]. Briefly, PCR was performed in 10 µL reactions containing 2 µL DNA template (around 50 ng), 1 µL of 10X TB buffer (containing 0.2 µL dd water; 0.5 µL of 1 M Tris-HCl, pH 8.4; 0.2 µL of 5 M KCl; 0.1 µL of 15 mM MgCl 2 ; 0.00000001 g gelatin), 0.5 µL of 1 mMdNTPs, 0.5 µL each of 10 µM forward and reverse primers, 5.3 µL nano pure ddH 2 O, and 0.2 µLTaq DNA polymerase (5 U/µL) using a thermal cycler. After initial denaturation for 2 min at 94 • C, each cycle comprised 30 s of denaturation at 94 • C, 30 s of annealing at 55 • C, and 30 s of extension at 72 • C with a final extension for 5 min at 72 • C at the end of 32 cycles. After completion of the reaction, the PCR product was preserved at a temperature of 4 • C. The amplified PCR products were separated in 6% polyacrylamide gel at 100 V for 1.5 to 2.5 h against 1Kb DNA marker. DNA bands were visualized with the exposure to ultraviolet light in a gel documentation system and saved in jpeg format. Allele scoring was performed considering the relative position of the bands in the gel images compared to the position of parental bands. Summary statistics, including the number of alleles per locus, major allele frequency, and polymorphic information content (PIC) values, were determined using Power Marker version 3.25 [41] following the formula proposed by Anderson et al. [42].
where n is the number of marker alleles for marker i and Pij is the frequency of the jth allele for marker i. Spearman's correlation coefficient was analyzed using the cor.test() function in R studio. The formula for calculating the Spearman correlation is as follows: where r s = Spearman correlation coefficient; di = the difference in the ranks given to the two variables values for each item of the data; n = total number of observations. The genetic distance was calculated using Nei distance [43]. The similarity matrix was calculated with the Simqual sub-program using the Shannon coefficient [44] and subjected to cluster analysis by the unweighted pair group method for the arithmetic mean (UPGMA), and a dendrogram was generated using the program NTSYS-pc [45].

Correlation among LD Score, Survivability, and Recovery Rate
The Spearman correlation among the cold-responsive traits measured in this study showed that the LD score, survivability, and recovery rate were significantly correlated to each other ( Figure 2). The highest correlation coefficient was observed between LD score and recovery rate, and the lowest correlation was between LD score and survival rate. However, both were in a negative or reverse direction. The survival rate and recovery rate were positively correlated. The principal component analysis also showed that LD score was negatively as-  The highest correlation coefficient was observed between LD score and recovery rate, and the lowest correlation was between LD score and survival rate. However, both were in a negative or reverse direction. The survival rate and recovery rate were positively correlated. The principal component analysis also showed that LD score was negatively associated with survivability and recovery rate (Figure 3).

Figure 2.
Spearman correlation between LD score, survivability, and recovery rate. The upper pan describes the correction coefficients and their significance level, the lower pan shows the directions of correlation, and the diagonal pan shows the distributions of the trait value across the genotypes.
The highest correlation coefficient was observed between LD score and recovery rate, and the lowest correlation was between LD score and survival rate. However, both were in a negative or reverse direction. The survival rate and recovery rate were positively correlated. The principal component analysis also showed that LD score was negatively associated with survivability and recovery rate (Figure 3).

Evaluation of NILs for Agronomic Performance
The NILs and the check varieties showed distinct genetic variation, particularly in grain yield per plant, the number of fertile grains per panicle, and % spikelet fertility. Days to flowering, plant height, and panicle length were also variable. Fertile grain per panicle had the widest range, followed by spikelet fertility. IR90688-92-1-1-1-1-1 had the highest value and IR90688-52-1-1-1-1-1 had the lowest value of FGPP. The ranges of panicle length, spikelet fertility, and thousand-grain weight were smaller compared to other characteristics. The maximum coefficient of variation was found in the case of fertile grain (12.7%) followed by grain yield per plant (11.43%).
The mean values of the NILs varied significantly from the check varieties for the characters measured in this study ( Table 2). The check varieties also differed from each other. There were many NILs that exceeded the mean values of the four high-yielding check varieties. Eleven NILs showed a significantly higher yield than the susceptible check varieties BR1 and BRRI dhan28. Out of these 11 NILs, IR90688-92-1-1-1-1-1, IR90688-125-1-1-1-1-1, IR90688-104-1-1-1-1-1, IR90688-124-1-1-1-1-P2, IR90688-15-1-1-1-1-1, and IR90688-27-1-1-1-1-1 had a significantly shorter growth duration than BR1, while only IR90688-92-1-1-1-1-1 and IR90688-125-1-1-1-1-1 were on par with BRRI dhan28 in terms of days to maturity. The NILs in this study were at an acceptable range in plant height except for one line ( Table 1). The genotypes showing significantly higher yield than BR18 and BRRI dhan28 had an intermediate type plant height. Not only that, these lines had an acceptable range of ineffective tillers per hill and filled grain per panicle (82-160), but also had a good seed setting rate under natural field conditions. In addition, they had almost a similar amount of thousand-grain weight to that of BRRI dhan28 and BRRI dhan29, which are the most widely grown varieties in Bangladesh. Principal component analysis of the agronomic traits with the cold-related seedling traits showed that grain yield was highly associated with PH and FG than other traits measured in this study.

Genetic Analysis of NILs of BRRI dhan29 Using Microsatellite Markers
Out of 88 SSRs used in this study, 58 markers were found polymorphic, with an average polymorphism rate of 36.97%. Among the 58 SSRs, 11 markers were present on chromosome 1, 6 markers on each of chromosomes 2 and 3, 4 markers on each of chromosomes 4, 9, and 12, 2 markers on each of chromosomes 5 and 6, 5 markers on chromosome 7, 3 markers on each of chromosome 8 and 10, and 8 markers on chromosome 11. A total of 186 alleles of the polymorphic SSRs were detected in the 42 genotypes, including 36 NILs (Table 3).  The number of alleles generated by each marker varied from 2 to 5 alleles and with an average of 3.21 alleles per locus. The highest number of alleles (5) was detected with RM12769 and the lowest number (2) was detected in RM10800, RM16686, RM24087, RM26063, RM26501, RM512, and RM7558. The overall size of the PCR products of the tested markers ranged from 72 bp (RM413) to 551 bp (RM13155). The molecular size difference between the smallest and the largest allele for a given locus varied widely from 3 to 112 bp (Table 3). There was a considerable range in allele frequency (42.86-83.33%). On average, 57.31% of the genotypes shared a common allele.
The PIC value, which indicates the degrees of allelic diversity [46] and frequency among the varieties, was calculated for all of the markers. The PIC values for the microsatellite loci varied from 0.27 to 0.54 with a mean of 0.42. Among the markers used in the study, RM553 on chromosome 9 had the highest PIC value (0.54). RM581 and RM11570 on chromosome 1, RM3703 on chromosome 2, RM14795 on chromosome 3, and RM264 on chromosome 8 showed PIC values greater than 0.50.
The dendrogram constructed based on Shannon's similarity index following the UPGMA method (Figure 4) showed that all of the genotypes were constellated into two distinct groups at a similarity index of 0.60. The first group housed 40 genotypes, while the second contained only 2 genotypes including an NIL and a cold-tolerant check variety, Hbj.B.VI. However, at the 0.67 similarity index, the genotypes were grouped into four major clusters, although the genotypes within the clusters were further subdivided into many subclusters based on the genetic similarity between them.
The dendrogram constructed based on Shannon's similarity index following the UP-GMA method (Figure 4) showed that all of the genotypes were constellated into two distinct groups at a similarity index of 0.60. The first group housed 40 genotypes, while the second contained only 2 genotypes including an NIL and a cold-tolerant check variety, Hbj.B.VI. However, at the 0.67 similarity index, the genotypes were grouped into four major clusters, although the genotypes within the clusters were further subdivided into many subclusters based on the genetic similarity between them.
Plant height in rice is an important trait for the better adoption of a variety. Earlier study revealed that the taller plants usually produce a low yield [53]; through another study reported that short-statured plants also give a low yield some extends due to the production of low biomass [54]. The NILs in this study were at an acceptable range of plant height, except for one line ( Table 1). The genotypes showing significantly higher yield than BR18 and BRRI dhan28 had intermediate type plant height. Not only did these lines have an acceptable range of effective tillers per hill and filled grain per panicle (82-160), but they also had a good seed setting rate under natural field conditions. In addition, they had an almost similar amount of thousand-grain weight to that of BRRI dhan28 and BRRI dhan29, which are the most widely grown varieties in Bangladesh. The principal component analysis of the agronomic traits with the cold-related seedling traits showed that grain yield was more highly associated with PH and FG than other traits measured in this study.
In addition to measuring yield potential, genetic diversity in the parental lines is very important for the genetic improvement of rice through breeding. Strong genetic diversity means diverse morphological traits and potentially valuable genetic information. Rice varieties/lines with high levels of genetic variation are essential resources for broadening the genetic base of the breeding germplasm, and therefore laying a good foundation for rice breeding [55]. Hence, the assessment of genetic diversity becomes important in establishing relationships among different cultivars. SSR markers are widely used in the analysis of genetic diversity among the germplasms [56].
Usually, NILs are developed through multiple backcrossing, meaning they are very close to each other as they share the majority of the recipient parental genome. In this study, we used BC 3 F 5 NILs developed from a cross between a Korean japonica variety Jinbubyeo and a Bangladeshi high-yielding rice variety BRRI dhan29, but they showed considerable variations in different agronomic traits. We analysed the NILs along with their parents, cold-tolerant and susceptible varieties, using 58 polymorphic SSR markers randomly distributed over 12 chromosomes to estimate the genetic similarity/dissimilarity between them. A total of 186 alleles of the polymorphic SSRs were detected in the genetic analysis of the 42 genotypes, including 36 NILs (Table 3). The number of alleles generated by each marker varied from 2 to 5 alleles and with an average of 3.21 alleles per locus. The highest number of alleles (five) was detected in the locus RM12769 and the lowest number of alleles (two) was detected in seven SSR loci, namely RM10800, RM16686, RM24087, RM26063, RM26501, RM512, and RM7558. The overall size of the PCR products of the tested markers ranged from 72 bp (RM413) to 551 bp (RM13155). The molecular size difference between the smallest and the largest allele for a given locus varied widely, from 3 to 112 bp ( Table 3).
The PIC value, which indicates the degrees of allelic diversity [46], varied from 0.27 to 0.54 with a mean of 0.42. The higher the PIC value, the greater the diversity within the germplasm; contrarily, a lower PIC value indicates that less divergence between the genotypes and they are closely related [40,57] in terms of the specific marker loci. Among the markers used in this study, RM553 on chromosome 9 had the highest PIC value (0.54). Two SSRs (RM581, RM11570) on chromosome 1, one SSR (RM3703) on chromosome 2, one SSR (RM14795) on chromosome 3, and one SSR (RM264) on chromosome 8 showed PIC values greater than 0.50, indicating that these are highly informative markers. An informative marker can easily and clearly discriminate the germplasm based on its band intensity and position relative to others [58]. Thus, highly informative markers are considered to be the best markers for diversity analysis. The six informative SSR markers identified from this study could be useful for future genetic studies of rice, particularly for marker-assisted breeding and QTL identification.
The UPGMA dendrogram constructed based on Shannon's similarity index (Figure 4) grouped all of the genotypes into two distinct clusters at a similarity index of 0.60. However, at the 0.67 similarity index, the genotypes were grouped into four major clusters, although the genotypes within the clusters were further subdivided into many subclusters based on the genetic similarity between them. Cluster 1, Cluster 2, Cluster 3, and Cluster 4 housed 25, 11, 4, and 2 genotypes, respectively. Housing the maximum number of genotypes, Cluster 1 was the most diverged cluster. The genotypes within Cluster 1 were further divided into many subgroups until 95% genetic similarity was reached. At 95% genetic similarity, IR90688-54-1-1-1-1-1 (G12) and IR90688-56-1-1-1-1-1 (G13) were grouped together. The genetic similarity among the genotypes in Cluster 1 ranged from 67.5% to 95%. The majority of the cold-tolerant and high-yielding NILs were grouped into Cluster 1. Contrary, Cluster 2 contained 11 genotypes, which were mostly poor yielders except for BR1 (G37), BR18 (G38), and BRRI dhan29 (G40). In this cluster, both cold-tolerant and susceptible varieties remained together, but they entered into different subgroups at a 0.77 similarity index. The most interesting thing is that Hbj. B.VI (G42), a cold-tolerant local variety for both seedling and reproductive stage [26], and IR90688-124-1-1-1-1-P2 (G35), which showed moderate tolerance to cold stress at the seedling stage (Table 1) and high yield (Table 2), were grouped into Cluster 4 with approximately 68% genetic similarity. Therefore, the crosses between IR90688-124-1-1-1-1-P2 and Hbj.B.VI could produce high yielding and cold-tolerant progenies. Not only these combinations, but also the cross combination between the lines of Cluster 4 with the high-yielding and cold-tolerant lines in Cluster 1, would produce cold-tolerant and high-yielding progeny.