*Article* **Phenotypic and Molecular Characterization of Rice Genotypes' Tolerance to Cold Stress at the Seedling Stage**

**Nasira Akter <sup>1</sup> , Partha Sarathi Biswas <sup>2</sup> , Md. Abu Syed 2,\* , Nasrin Akter Ivy <sup>3</sup> , Amnah Mohammed Alsuhaibani <sup>4</sup> , Ahmed Gaber <sup>5</sup> and Akbar Hossain 6,\***


**Abstract:** 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 artificially 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 significant 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 significantly 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.

**Keywords:** low-temperature stress; leaf discolouration; seedling stage; SSR markers; near-isogenic line; rice (*Oryza sativa* L.)

## **1. Introduction**

Compared to other cereal crops, rice is more sensitive to low-temperature stress (LTS) as it has originated from tropical regions [1–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–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

**Citation:** Akter, N.; Biswas, P.S.; Syed, M.A.; Ivy, N.A.; Alsuhaibani, A.M.; Gaber, A.; Hossain, A. Phenotypic and Molecular Characterization of Rice Genotypes' Tolerance to Cold Stress at the Seedling Stage. *Sustainability* **2022**, *14*, 4871. https://doi.org/10.3390/ su14094871

Academic Editor: Balázs Varga

Received: 20 March 2022 Accepted: 17 April 2022 Published: 19 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

early January [11]. Farmers are bound to replant the seedlings and, as a result, the cost of cultivation is increased. Sometimes, farmers delay the planting of the boro rice crops until the ambient temperature increases. Thus, they delay the next crop, which in turn results in the low total productivity of the farm. Moreover, it increases the risk of crop loss due to flash floods during the ripening stage in some low-lying areas of the north-eastern haors [4]. On the other hand, 100% spikelet sterility is observed in the short-duration rice varieties due to LTS at the booting stage of the rice crop during early to mid-February [12]. This problem could be resolved easily by developing cold-tolerant rice varieties. The success of a rice-breeding program depends on the presence of genetic variability among the parental lines used in it.

Genetic similarity/diversity can be assessed by using plant morphology, physiology, isozymes, storage protein profiles, DNA markers, etc. [13]. Among the different DNA markers, SSR markers have been extensively used as a powerful tool in variety protection [14], molecular diversity studies [15], quantitative trait loci (QTL) analysis, pedigree analysis, and marker-assisted breeding [16]. Morphological differences among near-isogenic lines (NILs) are very narrow [17]. The characterization of such a population is unrealistic if only based on morphological traits. Genetic analysis using molecular markers could reveal the difference between them at the DNA level. The molecular analysis of genes or QTLs underlying cold tolerance is the best approach, along with phenotypic cold tolerance screening. There are many reports available on QTL's underlying cold tolerance in the literature [18–26]; however, most of these QTL were detected from *japonica* cultivars.

A wide range of variations in cold tolerance exists among rice germplasms [19,27] and, overall, *Japonica* rice germplasms are more tolerant to cold stress than the *indica* rice germplasm [28–30]. To enhance cold tolerance in the *indica* rice, the commonly grown rice varieties in south and southeast Asian countries, it is imperative to use germplasms from temperate rice-growing countries as donor parents in the breeding program [18]. Crosses between *indica* and *japonica* rice usually show spikelet sterility due to cross incompatibility. However, marker-assisted backcrossing can recover recurrent parental genomes successfully. NILs, thus developed, could be used directly as a cold-tolerant variety after necessary evaluation or as the cold-tolerant parent as a donor to introgress into an elite *indica* background. In this study, we evaluated a set of NILs derived from a cross between a cold-tolerant *japonica* rice variety, Jinbubyeo and BRRI dhan29, a mega rice variety [31] of the boro ecosystem, to select potential lines for direct varietal release or use in the breeding program as the donor parents for cold tolerance.

#### **2. Materials and Methods**

#### *2.1. Sources of Materials Used*

Thirty-six BC3F<sup>5</sup> 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.

#### *2.2. 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.

BR18 and BR1 were used as tolerant and susceptible check varieties, respectively, in each tray (Figure 1). BR18 and BR1 were used as tolerant and susceptible check varieties, respectively, in each tray (Figure 1).

**Figure 1.** A pictorial view of the cold-screening nursery by using a cold water tank in the coldscreening laboratory of the Bangladesh Rice Research Institute (BRRI). **Figure 1.** A pictorial view of the cold-screening nursery by using a cold water tank in the coldscreening laboratory of the Bangladesh Rice Research Institute (BRRI).

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: ̇ 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:

$$\text{Survivaltivity rate } (\%) = \frac{\dot{\text{No. of Green plants}} \times 100}{\text{No. of plants treated}} \tag{1}$$

No. of Green plants × 100

(1)

$$\text{Recovery rate} \left( \% \right) = \frac{\text{No. of recovered plants} \times 100}{\text{No. of plants treated}} \tag{2}$$

mented RCB design [33]. At 35 days old, seedlings were transplanted with a single seed-

0.92 package in R 3.6.2. *2.3. Evaluation of Near Isogenic Lines (NILs) for Agronomic Performance* All 36 NILs along with four check varieties were grown in the field following aug-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.

#### ling per hill at a spacing of 20 cm × 20 cm in 5.4 m<sup>2</sup> plots. The standard crop management *2.3. Evaluation of Near Isogenic Lines (NILs) for Agronomic Performance*

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 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 <sup>×</sup> 20 cm in 5.4 m<sup>2</sup> 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:

$$R\mathbf{j} = B\mathbf{j} - M\tag{3}$$

where, *Rj* = effect of jth block, *Bj* = mean of all check-in the *j*th 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.

#### *2.4. 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 MgCl2; 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 ddH2O, 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].

$$\text{PICi} = \sum\_{j=1}^{n} (Pij)^2 \tag{4}$$

where *n* is the number of marker alleles for marker *i* and *Pij* is the frequency of the *j*th 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:

$$\mathbf{r\_s} = 1 - \frac{6\sum di\mathbf{2}}{n(n\mathbf{2} - 1)}\tag{5}$$

where r<sup>s</sup> = 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].

#### **3. Results**

#### *3.1. Evaluation of NILs for Cold Tolerance at Seedling Stage*

The analysis of variance for LD score, survivability, and recovery rate showed highly significant variation (*p* ≤ 0.01) among the genotypes (Table 1).

**Table 1.** Variance and mean performance of 39 genotypes in LD score, % survivability, and % recovery under artificial cold stress.


\* and \*\* indicate significant difference from BR18 at 5% and 1% level of probability, respectively.

#### 3.1.1. LD Score

LD scores among the genotypes varied from 2.2 to 8.2, with an average of 3.49 and CV of 15.17 (Table 1). The highest LD value was observed with the cold-susceptible check variety BR1 (8.2), which was statistically similar to those of IR90688-120-1-1-1-1-1, IR90688-124-1-1-1-1-P2, and BRRI dhan28 at the 5% level of significance. On the other hand, IR90688-103-1-1-1-1-1, IR90688-105-1-1-1-1-1, IR90688-5-1-1-1-1-1, IR90688-105-1-1-1-1-1, IR90688-118-1-1-1-1-1, IR90688-81-1-1-1-1-1, IR90688-20-1-1-1-1-1, and IR90688-74-1-1-1-1- 1 showed significantly lower LD values than that of BR18, which showed a moderate level of cold tolerance (LD score of 5.0). In addition, IR90688-62-1-1-1-1-1, IR90688-42-1-1-1-1-p2, IR90688-96-1-1-1-1-1, IR90688-54-1-1-1-1-1, and BRRI dhan29 showed LD values ranging from 3.9 to 4.5, indicating a moderate level of cold tolerance at the seedling stage.

## 3.1.2. Survivability Rate

The survivability rate varied widely among the genotypes, ranging from 38.70% to 100% with 19.20% CV. Out of 39 genotypes, 31 showed more than 60% survivability, while susceptible variety BR1 showed only 38.7% survivability (Table 1). IR90688-81-1-1-1-1- 1 showed the highest level of survivability (100%), followed by IR90688-103-1-1-1-1-1 (98.3%), IR90688-74-1-1-1-1-1 (95.0%), and IR90688-94-1-1-1-1-1 (90.9%), which were, in fact, not significantly different from the cold-tolerant check variety BR18 (76.9%), but were significantly different from the susceptible check variety BR1.

#### 3.1.3. Recovery Rate

The recovery rate also showed a wide range of variations among the genotypes, starting from 31.5% to 92.6% with an average value of 69.07 (Table 1). Although the tolerant and susceptible check varieties showed no significant difference in recovery rate, 23 NILs showed significantly higher recovery than BR1, and 4 (IR90688-74-1-1-1-1-1, IR90688-103- 1-1-1-1-1, IR90688-118-1-1-1-1-1, and IR90688-81-1-1-1-1-1) showed significantly higher recovery than even the cold-tolerant check variety BR18.

#### 3.1.4. 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). *Sustainability* **2022**, *14*, x FOR PEER REVIEW 7 of 17

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

related. The principal component analysis also showed that LD score was negatively as-

**Figure 3.** PCA plot showing an association between agronomic and cold-related seedling traits. The traits representing the same quadrant and in the same direction are positively associated, and the

sociated with survivability and recovery rate (Figure 3).

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

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 7 of 17

**Figure 3.** PCA plot showing an association between agronomic and cold-related seedling traits. The traits representing the same quadrant and in the same direction are positively associated, and the **Figure 3.** PCA plot showing an association between agronomic and cold-related seedling traits. The traits representing the same quadrant and in the same direction are positively associated, and the traits in the opposite direction are negatively associated. The angular distance between traits denotes the strength of association.

#### *3.2. 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.


**Table 2.** Adjusted mean of yield and yield-contributing traits of Jinbubyeo-BRRI dhan29 NILs.

GD = growth duration, PH = plant height, PL = panicle length, ET = effective tiller, FG = fertile grain, SF = spikelet fertility, TGW = thousand-grain weight, Ck = check variety, CV = coefficient of variation.

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.

#### *3.3. 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).

**Table 3.** Allelic variation and polymorphism information content values of different rice cultivars for 58 SSR polymorphic markers.



**Table 3.** *Cont.*

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 genetic similarity within the individual cluster/group ranges from 0.60 to 0.96. 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. Cluster 1 was subdivided into 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 17, 18, 20, 23, 24, and 25 subgroups at a genetic similarity of 67.5, 68, 70, 71, 71.5, 72.75, 73, 73.5, 74.5, 75, 76, 77, 78, 78.5, 79.5, 85, 91, and 95%, respectively. Contrarily, Cluster 2 contained 11 genotypes. The genetic similarity within Cluster 2 varied from 68 to 92%. IR90688-42-1-1-1-1-P1 (G8), IR90688-42-1-1-1-1-P2 (G9), IR90688-73-1-1-1-P1 (G16),

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 11 of 17

average, 57.31% of the genotypes shared a common allele.

chromosome 8 showed PIC values greater than 0.50.

and IR90688-82-1-1-1-1-1 (G20) were housed in Cluster 3. In Cluster 4, two genotypes were grouped at 67.5% genetic similarity. 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.

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

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

The PIC value, which indicates the degrees of allelic diversity [46] and frequency

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,

**Figure 4.** Dendrogram showing the genetic relationships among 42 rice cultivars using UPGMA cluster analysis of Shanondice genetic similarity coefficients generated from 58 SSR markers. G01 = IR90688-5-1-1-1-1-1; G02 = IR90688-13-1-1-1-1-P1; G03 = IR90688-15-1-1-1-1-1; G04 = IR90688-19-1-1- 1-1-1; G05 = IR90688-20-1-1-1-1-1; G06 = IR90688-27-1-1-1-1-1; G07 = IR90688-30-1-1-1-1-1; G08 = IR90688-42-1-1-1-1-P1; G09 = IR90688-42-1-1-1-1-P2; G10 = IR90688-43-1-1-1-1-P1; G11 = IR90688-52- 1-1-1-1-1; G12 = IR90688-54-1-1-1-1-1; G13 = IR90688-56-1-1-1-1-1; G14 = IR90688-62-1-1-1-1-1; G15 = IR90688-64-1-1-1-1-1; G16 = IR90688-73-1-1-1-1-P1; G17 = IR90688-74-1-1-1-1-1; G18 = IR90688-77-1- 1-1-1-1; G19 = IR90688-81-1-1-1-1-1; G20 = IR90688-82-1-1-1-1-1; G21 = IR90688-91-1-1-1-1-1; G22 = IR90688-92-1-1-1-1-1; G23 = IR90688-94-1-1-1-1-1; G24 = IR90688-95-1-1-1-1-1; G25 = IR90688-96-1-1- 1-1-1; G26 = IR90688-103-1-1-1-1-1; G27 = IR90688-104-1-1-1-1-1; G28 = IR90688-105-1-1-1-1-1; G29 = IR90688-106-1-1-1-1-1; G30 = IR90688-108-1-1-1-1-1; G31 = IR90688-109-1-1-1-1-1; G32 = IR90688-114- 1-1-1-1-1; G33 = IR90688-118-1-1-1-1-1; G34 = IR90688-120-1-1-1-1-1; G35 = IR90688-124-1-1-1-1-P2; G36 = IR90688-125-1-1-1-1-1; G37 = BR1; G38 = BR18; G39 = BRRI dhan28; G40 = BRRI dhan29; G41 = Jinbubyeo; G42 = Hbj.B.VI. **Figure 4.** Dendrogram showing the genetic relationships among 42 rice cultivars using UP-GMA cluster analysis of Shanondice genetic similarity coefficients generated from 58 SSR markers. G01 = IR90688-5-1-1-1-1-1; G02 = IR90688-13-1-1-1-1-P1; G03 = IR90688-15-1-1-1-1-1; G04 = IR90688-19-1-1-1-1-1; G05 = IR90688-20-1-1-1-1-1; G06 = IR90688-27-1-1-1-1-1; G07 = IR90688- 30-1-1-1-1-1; G08 = IR90688-42-1-1-1-1-P1; G09 = IR90688-42-1-1-1-1-P2; G10 = IR90688-43-1-1-1-1-P1; G11 = IR90688-52-1-1-1-1-1; G12 = IR90688-54-1-1-1-1-1; G13 = IR90688-56-1-1-1-1-1; G14 = IR90688- 62-1-1-1-1-1; G15 = IR90688-64-1-1-1-1-1; G16 = IR90688-73-1-1-1-1-P1; G17 = IR90688-74-1-1-1-1-1; G18 = IR90688-77-1-1-1-1-1; G19 = IR90688-81-1-1-1-1-1; G20 = IR90688-82-1-1-1-1-1; G21 = IR90688- 91-1-1-1-1-1; G22 = IR90688-92-1-1-1-1-1; G23 = IR90688-94-1-1-1-1-1; G24 = IR90688-95-1-1-1-1-1; G25 = IR90688-96-1-1-1-1-1; G26 = IR90688-103-1-1-1-1-1; G27 = IR90688-104-1-1-1-1-1; G28 = IR90688-105- 1-1-1-1-1; G29 = IR90688-106-1-1-1-1-1; G30 = IR90688-108-1-1-1-1-1; G31 = IR90688-109-1-1-1-1-1; G32 = IR90688-114-1-1-1-1-1; G33 = IR90688-118-1-1-1-1-1; G34 = IR90688-120-1-1-1-1-1; G35 = IR90688- 124-1-1-1-1-P2; G36 = IR90688-125-1-1-1-1-1; G37 = BR1; G38 = BR18; G39 = BRRI dhan28; G40 = BRRI dhan29; G41 = Jinbubyeo; G42 = Hbj.B.VI.

#### **4. Discussion**

Cold tolerance at the seedling and reproductive stage is very important for growing rice in the winter season or in cold stress conditions. The early vegetative stages, particularly at the germination and crop establishment stage and also the booting stages of rice, are very sensitive to LTS (Figure 1). An ambient temperature below 10–12 ◦C significantly increases seedling mortality in rice [47], resulting in poor crop establishment, which eventually turns into low yield. By assigning LD scores [48] upon exposure to LTS for seven days in artificially induced low temperature [32] and measuring the survivability and recovery rate, the genotypes with seedling-stage cold tolerance can easily be identified. In this study, we evaluated 36 NILs with cold tolerance and susceptible check varieties for seedlingstage cold tolerance. The genotypes showed a wide variation in LD scores, ranging from 2.2 to 8.2 with an average of 3.49 (Table 1). Suh et al. [49] also reported varying LD scores (3 to 8 with a mean value of 5) among 23 elite rice cultivars grown in a cold waterirrigated plot. In this study, the highest LD value was found with the cold-susceptible check variety BR1 (8.2). IR90688-120-1-1-1-1-1, IR90688-124-1-1-1-1-P2, and BRRI dhan28 showed

no significant difference from BR1 in LD score at *p* < 0.05, indicating them to be coldsusceptible. Contrarily, IR90688-103-1-1-1-1-1, IR90688-105-1-1-1-1-1, IR90688-5-1-1-1-1-1, IR90688-105-1-1-1-1-1, IR90688-118-1-1-1-1-1, IR90688-81-1-1-1-1-1, IR90688-20-1-1-1-1-1, and IR90688-74-1-1-1-1-1 showed significantly lower LD values than the cold-tolerant check variety BR18. These results clearly revealed the capability of tolerance of the NILs to LTS. Additionally, IR90688-62-1-1-1-1-1, IR90688-42-1-1-1-1-p2, IR90688-96-1-1-1-1-1, IR90688- 54-1-1-1-1-1, and BRRI dhan29 showed LD values ranging from 3.9 to 4.5, indicating a moderate level of cold tolerance at the seedling stage. Khatun et al. [32] reported a higher average LD score for BR1, moderate for BRRI dhan29, and the lowest LD score for BR18. Kundu [50] also reported similar findings from a cold-screening trial of a set of cold-tolerant *japonica* and *indica* donors, along with cold-susceptible cultivars.

The seedling survivability and recovery rate are very important for better crop establishment under LTS. In this study, the survivability and recovery rate varied widely among the genotypes, ranging from 38.70% to 100% and 31.5% to 92.6%, respectively. The susceptible variety BR1 showed only 38.7% survivability and 31.5% recovery rate (Table 1). IR90688-81-1-1-1-1-1 showed the highest level of survivability (100%), followed by IR90688- 103-1-1-1-1-1 (98.3%), IR90688-74-1-1-1-1-1 (95.0%), and IR90688-94-1-1-1-1-1 (90.9%), which were, in fact, not significantly different from the cold-tolerant check variety BR18 (76.9%), but were significantly different from the susceptible check variety BR1. A total of 23 NILs showed a significantly higher recovery rate than BR1, with 4 of them (IR90688-74-1-1-1-1-1, IR90688-103-1-1-1-1-1, IR90688-118-1-1-1-1-1, and IR90688-81-1-1-1-1-1) showing significantly higher recovery than even the cold-tolerant check variety BR18. In a similar study with 10 landraces and 5 modern varieties as check genotypes, Shoroardi et al. [51] showed a higher reduction in seedling survival rate in BRRI dhan28, LR 212- NariaBochi, and LR 25-KathiGoccha and a lower reduction in BRRI dhan36, BRRI dhan29, LR 54- Tor Balam, and LR 114-Badiabona. The magnitude of the recovery rate was a little lower than the survivability rate, but they were positively correlated. The recovery rate, in fact, pronounced the real tolerance to LTS at the seedling stage.

For direct use in variety release or as a parental line in the breeding program, a breeding line should have higher yield potential in addition to its particular trait of interest. Therefore, we evaluated the NILs for yield and other agronomic performances under natural field conditions. Although the NILs were developed in the background of BRRI dhan29, they possessed distinct genetic variations, particularly in grain yield per plant, number of fertile grains per panicle, and % spikelet fertility [52]. Days to flowering, plant height, and panicle length were also variable among the NILs. The mean values of the NILs varied significantly from the check varieties for all of the characteristics measured in this study (Table 2). Eleven NILs showed significantly higher yields than the susceptible check varieties BR1 and BRRI dhan28. Out of these 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.

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 BC3F<sup>5</sup> 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.

#### **5. Conclusions**

Highly significant variation among the genotypes was observed for three cold-related traits at the seedling stage, namely, leaf discolouration score, survivability, and recovery rate. Nine NILs showed strong cold tolerance in at least one of three cold-related traits. 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 to cold stress with all three traits, while IR90688-118-1-1-1-1-1 showed cold tolerance with LD score and recovery rate. In a field experiment, 11 NILs showed significantly higher yields 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 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. In addition, these lines had grain quality similar to BRRI dhan28 and BBRI dhan29. In the genetic analysis with 58 polymorphic SSRs, a total of 186 alleles were detected and the number of alleles per marker varied from 2 to 5, with an average of 3.21 alleles per locus. Six SSRs showed PIC values greater than 0.50, indicating that these are the best markers in discriminating the NILs derived from the cross between Jinbubyeo and BRRI dhan29. The Shannon similarity index-based dendrogram grouped the NILs at 67% genetic similarity into four major distinct clusters, in which 25 genotypes, including the high-yielding NILs and high-yielding variety BRRI dhan28, constellated together. One high-yielding NIL, IR90688-124-1-1-1-1-P2, with a moderate level of cold tolerance was grouped; however, genotype Hbj. Boro VI has been found to be a cold-tolerant local variety for both the seedling and reproductive stages, as it is a member of the two-genotype cluster. Therefore, the cold-tolerant and high-yielding NILs identified in this study could be a potential genetic resource to develop high-yielding, cold-tolerant rice varieties. The QTL mapping from the local variety Hbj. Boro VI may be considered to generate marker information, aiming to introgress into the locally adapted existing rice varieties for the sustainability of rice production during extremely cold stress conditions.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/su14094871/s1, Table S1: List of 88 SSR markers used for diversity analysis of rice cultivars.

**Author Contributions:** Conceptualization: N.A., N.A.I., M.A.S. and P.S.B.; methodology: N.A., N.A.I., M.A.S. and P.S.B.; formal analysis: N.A., N.A.I., M.A.S. and P.S.B.; data curation: N.A., M.A.S. and A.H.; statistical expertise: N.A., M.A.S. and A.H.; writing—original draft preparation: N.A., N.A.I., M.A.S. and P.S.B.; writing—review and editing: A.G., A.M.A. and A.H.; visualization: N.A., N.A.I., M.A.S. and P.S.B.; supervision: P.S.B.; funding acquisition: A.M.A., A.G. and A.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was financially supported by three institutes, namely, the Bangladesh Agricultural Research Institute, Bangabandhu Sheikh Mujibur Rahman Agricultural University, and the Bangladesh Rice Research Institute, Gazipur, Bangladesh. The research program also partially supported Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2022R65), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are available upon request.

**Acknowledgments:** The authors wish to acknowledge the financial support provided by Bangladesh Agricultural Research Institute, Bangabandhu Sheikh Mujibur Rahman Agricultural University, and the Bangladesh Rice Research Institute, Gazipur, Bangladesh. The authors appreciated Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2022R65), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia, for providing financial support for the publication of this research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

