Progeny (262 RILs) between Sherpa and IRAT109 were produced (F₁s on two separate occasions) by the NSW Department of Primary Industries (DPI). The resulting F₁s were then self-pollinated to produce F₂s and subsequently, single seed descent method was carried out until F₆ generation to produce the RILs. Black plastic tube pots were filled with Lockyer prairie soil consisting of light, black clay [USDA Soil Taxonomy: Fluventic Hapludolls; (Isbell, 2016)] across all the experiments carried out in this study. Fertilizers were added at a rate of 3 kg/m³ of slow release Osmocote^® 3–4M (19-9-10 + 2MgO + TE, ICL Specialty Fertilizers, Australia) and 2 kg/m³ of slow release Osmocote^® 5–6M (15-9-12 + 2MgO + TE) which translated to 0.72 and 0.48 g of each fertilizer per pot (0.24 L), respectively. Each pot was encased with a similar sized pot to prevent root escape. Approximately five seeds were sown in each pot and were subsequently thinned to a density of one seedling per pot 18 days after seeding (DAS). Seedlings were grown in an aerobic condition until 11 DAS, after which pots were gradually flooded with water level raised by 3 cm until the seedlings were fully flooded (3 cm above the soil surface). Iron sulfate and insecticide were sprayed as needed.

Alongside the 262 RILs, 26 checks were added in the experiments including known LT tolerant (Lijiangheigu and RL11) and susceptible (Reiziq and RL206) checks and other varieties and advanced breeding lines. The LT screening protocol was carried out as outlined by Mitchell et al. (2016) with some modifications. The 288 lines were seeded in two sets: Set 1 (Experiment 1) and Set 2 (Experiment 2). Experiment 1 was a partially replicated design. Experiment 2 was planted 16 days after Experiment 1 and was designed with three full replications. Both experiments were grown in a warm room (30°C/21°C day/night) with natural light. When a single plant of a particular genotype reached heading (when the first spikelet emerged from the sheath) in Experiment 1, the corresponding plants in Experiment 2 were moved to a low temperature room (20°C/15°C day/night) for 14 days for LT treatment. The plants were then moved back to the warm room until maturity. At maturity, plants in Experiment 2 were phenotyped for the following traits: plant height (PH), spikelet per panicle (SPP), and spikelet sterility (SS). Total number of spikelets were counted from mainstem panicles, for which filled (FS), unfilled, and dead spikelets were identified and counted independently. SS was determined using the ratio between the sum of unfilled and dead spikelets over the total number of spikelets multiplied by 100. Days to heading (DTH) for Experiment 1 genotypes were also recorded and analyzed.

The full population was re-evaluated for LT tolerance similar to the above experiments conducted with some modifications to confirm consistency of identified genomic regions associated with LT tolerance traits. Set 1 (Experiment 3) and Set 2 (Experiment 4) were seeded 16 days apart with two replications in each experiment in a temperature and light controlled glasshouse. Both sets were grown in a warm room (28°C/22°C day/night) and light was controlled with installation of black-out curtain from 5:00 PM to 7:00 AM. Corresponding plants from Experiment 4 were moved to a low temperature room (21°C/15°C day/night) when a single replicate from Experiment 3 reached heading and subsequently, the same steps were carried out as described in Experiments 1 and 2. DTH and SS from Experiments 3 and 4, respectively were collected and used for QTL mapping.

From the results of Experiments 1 and 2, a set of differential RILs which represented the phenotypic extremes in terms of SS of the population were selected to determine the responses of SS and floral traits to LT in Experiment 5 and 6. The RIL selection was based on DTH from the full population evaluation (<70, 70–75, 75–80, and >80 days after seeding), with tolerant and sensitive RILs having SS < 32% and SS > 68%, respectively. For each of these maturity groups, the top five and bottom five genotypes in terms of SS were selected. A further nine genotypes were added based on their relatively low SS from Experiment 2 (<32%) and along with two promising RILs (with other useful AP traits) bringing the total number of RILs to 51. Additionally, nine check genotypes (parents Sherpa and IRAT109, known LT tolerant donors – Norin PL8, Lijiangheigu, and RL11, and elite genotypes – Reiziq, YRE16 = V71, Viand, and YUA15 = V037) were also included. Using the two-set system described above, these 60 genotypes each with three replications were grown in warm conditions throughout the growing period (30°C/21°C day/night) in Experiment 5, while in Experiment 6 they were exposed to low temperature conditions (20°C/15°C day/night) at YMS for 14 days.

Spikelet sterility, DTH, and PH of the genotypes from Experiments 5 and 6 were recorded. Main stems in all the plants tested were tagged and subsequent collection of spikelets for floral trait analyses were obtained from the panicle of the same stem. Floral traits including number of pollen on stigma (POS), number of dehisced anthers (NDA), and total number of pollen per spikelet (TP) using an impedance flow cytometer (Amphasys AG, Switzerland) were quantified. For the determination of TP, spikelets with mature anthers (about to open) were collected and processed. Briefly, six anthers per spikelet were dissected using a fine dissecting needle, after which 1.5 mL of AF6 buffer (Amphasys AG, Switzerland) was added and subsequently filtered using a 100 μm filter (Amphasys AG, Switzerland). The filtrate was then read using AmphaChip E00913, a 250 μm chip (Amphasys AG, Switzerland). Thresholds clearly separating pollen from debris was manually determined using the data from warm set-up and was subsequently used throughout the two experiments. TP was calculated by adding the upper left and upper right readings in the quadrant. NDA and POS were processed following methods described by Susanti et al. (2019) with some modifications. Shortly after floret opening, spikelets were collected from the middle of the main stem panicle. To determine POS, the stigma was dissected from the spikelets and the associated pollen grains were stained using a Lugol solution (Merck Pty. Ltd., Germany), and the number of fertile (engorged; stained black) were counted. The total NDA was counted out of the six anthers and along with POS, were visualized under a stereomicroscope (Olympus SZX10; Olympus Corporation, Japan). Using SS determined from Experiment 2 and this experiment, genotypes were then grouped into tolerant and sensitive RILs to facilitate statistical comparisons among the traits collected above.

Experiments were designed using the DiGGER package in R (Coombes, 2009). Statistical analysis was carried out on phenotypic traits collected from the full population (Experiments 1–4) using ASReml package in R (VSNi, United Kingdom). Best linear unbiased predictors (BLUPs) and heritabilities (Smith et al., 2006) of traits were computed using a model with genotype as a random effect. For the differential set experiment (Experiments 5 and 6), best linear unbiased estimators (BLUEs) were computed on a resolvable row-column design using a mixed linear model treating genotype as a fixed effect using the lme4 package in R (Bates et al., 2015). Testing for significance of difference between groups was also carried out using multcomp package in R (Hothorn et al., 2008). Significant genotypic variation among genotypes and groups tested were considered at p < 0.05, unless stated otherwise. Correlation between significant traits were determined through Pearson product-moment correlation coefficient using BLUPs and BLUEs computed in respective experiments.

Out of the 262 RILs, 252 were genotyped for single nucleotide polymorphism (SNP) using the Diversity Arrays Technology (DArT) genotyping-by-sequencing platform (DArTSeq). DArTSeq allowed for complexity reduction and sequencing of low copy sequences (corresponding predominantly to active genes) through the use of methylation sensitive restriction enzymes. These libraries were then sequenced using a next generation sequencing platform (Illumina HiSeq 2500) and the resulting sequences were aligned to Oryza sativa v7 reference genome to identify SNPs. DArTSeq resulted in a total of 6,508 SNP loci distributed across the rice genome and a total of 2,632 high quality polymorphic SNP loci were retained (Vinarao et al., 2021b) after filtering for call rates and markers following the expected allele frequency for RILs. Using the final genetic map and genotype data generated by Vinarao et al. (2021b) from the same population, a multiple QTL model (MQM) analysis (Broman et al., 2003) was carried out to identify genomic regions associated with SS and DTH collected from the full population. The logarithm of odds (LOD) threshold used to report QTL were based on the results of 1,000 permutations at a 5% significance level. Finally, using the bayesint function, Bayesian credible 95% confidence interval surrounding the identified QTL peak were determined in R/qtl (Broman et al., 2003).

High-quality genomic DNA from Sherpa was extracted using a modified cetyltrimethylammonium bromide extraction method (Doyle and Doyle, 1987). DNA quality was checked using 1% agarose gel electrophoresis and NanoDrop™ 8000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). Extracted DNA was sent to Macrogen Inc., South Korea where DNA libraries were prepared and sequenced. Library preparation was carried out using TruSeq Nano DNA 350 (Illumina, United States) and validated for concentration and size distribution using Agilent TapeStation D1000 Screen Tape (Agilent, United States). The resulting libraries were then sequenced with Illumina NovaSeq 6000 (Illumina, United States) paired-end 150bp reads. Raw reads (paired-end reads of 90 bp) of IRAT109 were retrieved from the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA: ERS470547) which was made publicly available through the 3K Rice Genotyping Project (3KRGP) (Mansueto et al., 2017).

Initial quality control analyses of the whole-genome re-sequencing (WGRS) were carried out using the FastQC tool (Andrews, 2010) and then Trimmomatic (Bolger et al., 2014) was used to remove Illumina adapter sequences, filter reads shorter than 36 bp, and to trim the start and end of the sequences having quality lower than 5 (Ignacio et al., 2021). The trimmed sequences were then aligned to the Oryza sativa MSU release 7 reference genome² using the Burrows-Wheeler Alignment maximal exact matches sequence aligner (Li, 2013), and were subsequently sorted and assigned respective groups. The resulting alignments were marked for PCR and optical duplicates and filtered for mapping quality (MAPQ) score greater than 20 (Li et al., 2009). Additionally, alignments with the following flag bits were retained: (1) read is paired and (2) read is mapped in a proper pair, and subsequent alignments were left realigned through the BamLeftAlign tool to homogenize positional distribution of insertions and deletions (indels). A bayesian genetic variant detection was carried out using FreeBayes to find SNPs, indels, multi-nucleotide polymorphisms and complex variants in a simple diploid calling manner with filtering and coverage: (1) MAPQ of 30, (2) minimum base quality of 20, and (3) minimum coverage of 15 to process a site (Garrison and Marth, 2012). Finally, variants were annotated and effect predicted using SnpEff tool (Cingolani et al., 2012). Variants were further filtered for polymorphism between parents using Trait Analysis by aSSociation, Evolution and Linkage software (Bradbury et al., 2007). Co-expression analysis for LOC_Os03g03100 was carried out using the online Rice Annotation Project Database (RAPDB) (Sakai et al., 2013) and gene ontology (GO) enrichment analysis was performed in RiceFREND (Sato et al., 2013). A built-in tool of the SNP-Seek database (Mansueto et al., 2017) was used to carry out haplotype analysis for LOC_Os03g03100 in the 3024 sequenced lines available through the 3KRGP. The Calinski criteria (Caliński and Harabasz, 1974) for automatic k-group determination was employed to cluster the lines into groups/haplotypes (Abbai et al., 2019). The All 3K SNP dataset containing 32 million SNPs was utilized.

In Experiment 1, 244 RILs established and were subsequently used as indicators for determining the time of transfer of the genotypes to LT glasshouse in Experiment 2. Highly significant genotypic variation (p < 0.0001) was observed for all the phenotypic traits collected in Experiment 2 with heritabilities ranging from 0.67 to 0.86 (Table 1). Ten genotypes were found to escape the LT treatment and as such a total of 234 RILs were included in the following further analysis. Lijiangheigu and RL11 had relatively low SS, with 46.7 and 21.6%, respectively. In terms of the parents, Sherpa produced a lower SS, compared to IRAT109, with 59.1 and 77.6%, respectively. Within the RILs, 8.6% (21 genotypes) had SS from 0-30%, and therefore were considered LT tolerant, while the majority of the RILs (129, 52.9%) had SS greater than 50% (Figure 1A). The final tolerant selected RILs (19) had an average sterility of 27.2%, while the susceptible RILs (20) had an average of 82.1%. In terms of PH, Lijiangheigu was the tallest (117.7 cm) among the checks and RILs tested. IRAT109 was also taller (77.8 cm) than Sherpa (65.2 cm). The tallest derivative of Sherpa/IRAT109 recorded PH of 98.4 cm while the shortest was 52.9 cm. There was also a significant genotypic variation in SPP. IRAT109 showed a slightly higher SPP compared with Sherpa, with 78.4 and 66.5, respectively. The RILs had SPP from as high as 124.3 to as low as 43.9. In terms of DTH collected in Experiment 1, highly significant genotypic variation (p < 0.001) was also observed, with a heritability of 0.94. Of the parents, Sherpa, headed earlier (65.8 DAS) than IRAT109 (75.6 DAS), while their progenies took 64.6–90.5 days to head (Figure 1A). Correlation analysis among the RILs revealed significant relationships among the traits collected. Earlier heading genotypes in Experiment 1 tended to have higher sterility in Experiment 2 (r = −0.34^**). Similarly, taller plants also tended to have lower SS (r = −0.24^**).

Highly significant genotypic variation (p < 0.0001) was also observed for DTH and SS collected in Experiment 3 and 4 with heritabilities of 0.73 and 0.81, respectively (Table 1). DTH in Experiment 3 showed an average of 74.9 days with the top 25 earliest RILs showing an average DTH of 66.9 while the top 25 latest RILs registered 81.6 days (Figure 1A). In terms of SS from Experiment 4, the average was 37.9% with the tolerant RILs producing an average SS of 27.1% while the sensitive RILs had 51.2%. SS between Experiments 2 and 4 were also significantly correlated (r = 0.42^**).

Highly significant genotypic variation (p < 0.001) was observed in the floral and agronomic traits collected in Experiment 5, except for NDA (Table 2). Heritabilities of the traits ranged from 0.67 (TP) to 0.96 (DTH). SS ranged from 0.7 to 25.7% with a mean of 5.8%. Comparing the parents, consistent with the full population evaluation, Sherpa headed earlier (66.4 days) than IRAT109 (71.6 days). In terms of the analysis of group effect, no significant difference was observed between the tolerant and sensitive RILs when grown under warm conditions—the groups had similar SS, were similar in days to heading, with similar floral trait results.

Highly significant genotypic variation (p < 0.001) was observed in the floral traits (NDA, POS, and TP) as well as the key agronomic traits [SS (Figure 1B), DTH, PH, and FS], when grown under LT in Experiment 6 with heritabilities ranging from 0.62 (POS) to 0.97 (DTH in LT experiment) (Table 2). RL11, Norin PL8, and Lijiangheigu had low SS with 9.1, 16.6, and 18.0%, respectively. Sherpa consistently produced lower SS (50.8%) compared with IRAT109 (69.9%). Their derived RILs showed an average of 51.0% SS, with genotypes as low as 2.3% to as high as 99.6%. Examining the difference between RIL genotype groups, tolerant RILs produced significantly lower SS with an average of 18.0% compared with the sensitive RILs (73.6%). Tolerant RILs showed significantly higher number of FS, more than double compared to the sensitive RILs. Analysis into the floral traits between these groups also indicated that tolerant RILs had significantly higher NDA, with an average 2.6 more dehisced anthers and an increased number of POS compared with sensitive RILs (92 vs 42, 119%). Furthermore, tolerant RILs had 72% more pollen per spikelet than the sensitive RILs. It is worth noting that there was no significant difference between DTH of tolerant and sensitive RIL groups owing to the stratified sampling strategy employed during genotype selection. In summary, after exposure to LT conditions at YMS, tolerant RILs had lower SS, produced higher NDA, POS, and TP compared to sensitive RILs. TP between Experiments 5 and 6 were related (r = 0.4^**), indicating that genotypes with high number of pollen per spikelet in warm conditions also tended to have higher TP in LT.

Of the 252 RILs with SNP genotype information, 241, 228, 214, and 229 RILs were used to map genomic regions associated with DTH in Experiment 1 and 3, and SS in Experiment 2 and 4, respectively. A total of ten RILs were determined to have escaped the LT treatment in Experiment 2 and were therefore excluded from the QTL mapping effort for SS. Using the corresponding SNP genotyping data and phenotypic data from the full population, linkage mapping was carried out and identified a total of nine genomic regions: three associated with SS while six were for DTH (Table 3, Figure 2, and Supplementary Figure 1). MQM analysis was carried out for SS with a LOD threshold of 3.06 and 3.12 and subsequently identified three and two genomic regions associated with lower SS for Experiments 2 and 4 respectively, explaining 19.5–22.4% of the total genotypic variation. Of the QTL identified, two (qYMCT3 and qYMCT8.1) were identified in both experiments while one (qYMCT4) was identified only in Experiment 2. qYMCT3 with an LOD score of 4.14–6.96 was located around 2.19 Mb physical position, with a confidence interval of 0.9–11.7 cM. Introgression of this QTL with the Sherpa (AA) allele was estimated to lower the SS of a genotype by 10.6–15.4% (Figure 3) and this region explained about 7.24–12.08% of the genotypic variation observed in the population. qYMCT8.1 with a LOD score of 3.07–3.26 was located around 18.45 Mb physical position of Chromosome 8. Introgression of the Sherpa allele from qYMCT8.1 was estimated to decrease SS of a genotype by 8.6–10.8% (Figure 3). The SS and SNP genotype composition of three selected RILs harboring the favorable (AA) allele of qYMCT3 and qYMCT8.1 (Supplementary Figure 2) were shown in Table 4 in contrast to the composition of one RIL harboring the unfavorable allele. Additionally, qYMCT4 was identified on Chromosome 4 (14.89 Mb) with a LOD score of 3.77 (Figure 2). In qYMCT4 the introgression of the IRAT109 (BB) allele decreased SS by about 9.6%.

As initial phenotypic analyses indicated the possible crosstalk between SS and DTH, MQM analysis was conducted and identified genomic regions associated with DTH in warm conditions, with LOD greater than the 3.09 and 3.08 computed threshold for Experiments 1 and 3, respectively. A total of six genomic regions located in five chromosomes (2, 3, 4, 5, and 12) were identified to influence DTH (Table 3, Figure 2, and Supplementary Figure 1). The final QTL model for Experiments 1 and 3 were able to explain 47.4–63.1% of the variation in DTH. In both experiments, qDTHW3.1 located between genetic distance of 0.9–7.6 cM was the major QTL identified in this population with a LOD score of 23.93–41.58 and alone was able to explain 32.7–44.7% of the genotypic variation. Interestingly, this QTL was shown to co-locate with qYMCT3 and the analysis also showed that introgression of the AA allele resulted in the genotypes to reach heading about 5.7–7.4 days (Figure 3) later compared to genotypes with BB allele. Of the five remaining QTL, three (qDTHW3.2, qDTHW4, and qDTHW5) were identified in both experiments while qDTHW2 and qDTHW12 were identified only in Experiment 1. The next major genomic region associated with DTH was found in chromosome 5 (qDTHW5) with a LOD score of 10.21–13.69 and an R² of 11–12.1%. Introgression of the Sherpa allele was shown to result in an earlier heading of the genotypes tested, up to about 3.8 days. qDTHW3.2 had a LOD score of 6.44–11.32 and explained 4.8–13.5% of the variation. Introgression of the Sherpa allele from qDTHW3.2 also resulted in an estimated 2.6–4.5 earlier DTH. The other remaining QTL (qDTHW2, qDTHW4 and qDTHW12) showed modest effects with LOD scores ranging from 3.75 to 9.54 and were able to explain between 2.7 and 11.2% of the variation. Similar to qDTHW5 and qDTHW3.2, introgression of the Sherpa allele from these QTL resulted in earlier heading with effects ranging from 2 3 days earlier than with IRAT109 allele. Significant interaction between qDTHW3.1 and qDTHW4 was found in Experiment 1 with a LOD score of 3.05 and explained 2.2% of the variation. Similarly, in Experiment 3, a significant epistatic interaction between qDTHW3.2 and qDTHW4 was also detected, with a LOD score of 3.72 and explaining 4.1% of the variation. Looking into the additive effects of the genomic regions identified, although qDTHW3.1 had the strongest effect (Sherpa allele cause delayed heading), the magnitude of the five QTL outweighed the effect of qDTHW3.1.

Using the MSU Rice Genome Annotation Project (Kawahara et al., 2013), a survey of previously identified genes associated with DTH along with the qDTHW3.1 confidence interval was carried out. It was found that two cloned genes associated with heading in rice: Ehd4 (LOC_Os03g02160, Chr03:717447-720837) and DTH3 (LOC_Os03g03100, Chr03: 1298070-1300384) were located within this region.

Analysis of the whole genome sequence data obtained from IRAT109 and Sherpa showed more than 109 and 69 million reads passed QC for the genotypes, respectively (Supplementary Table 1). Of these reads, a high percentage of the sequences were mapped with the reference Nipponbare genome at 95.8 and 98.2%, respectively for IRAT109 and Sherpa. After filtering for duplicate reads and retaining properly paired reads, the depth of coverage was calculated to be around 22X and 21X, for IRAT109 and Sherpa, respectively.

Variant calling carried out using Nipponbare reference sequence identified a total of 1,094,121 variants, of which 79.8, 7.65, 5.14, and 5.68%, were SNPs, multiple nucleotide polymorphisms (MNPs), insertions and deletions, respectively (Supplementary Table 2). Chromosome wise, chromosome 1 had the highest number of variants with 132,376 while chromosome 9 showed the lowest with 45,610. On average, the frequency of variant detection occurred every 342 number of bases. In terms of the number of variants by type and region, majority of the variants occurred in the intergenic region, followed by upstream and downstream gene variants, respectively. Finally, looking into the functional class of the variants, 55.6% of the variants were predicted to cause missense amino acid change in the protein product while 43.18 and 1.18% were predicted to cause silent and non-sense effects.

Using the variants identified, sequence variation between the parents in the two DTH genes identified to co-locate with the qDTHW3.1/qYMCT3 locus was investigated. It was found that no polymorphism occurred between Sherpa and IRAT109 in LOC_Os03g02160 while a total of 17 variants were identified in LOC_Os03g03100 gene (Figure 4 and Supplementary Table 3). Of the variants detected, 14 (82%) were SNPs while two were predicted to be complex variants and a single insertion variant. In terms of their effects, 10 were predicted to be upstream gene variants, one 5′ UTR variant, and the rest (six) were intron variants. Co-expression network analysis (hierarchy: 3, MR rank: 3) was carried out in RAPDB for LOC_Os03g03100 (RAPDB ID: Os03g0122600), and a total of 17 genes (Supplementary Table 4). Gene ontology (GO) enrichment analysis using these genes showed significant GO molecular functions (Supplementary Figure 3) related with calcium ion transporter activity (GO:0005388 and GO:0015085) as well as ATPase activity coupled with transmembrane movement of ions and substances (GO:0015662, GO:0042625, and GO:0042626). Haplotype analysis showed the presence of three haplotypes in LOC_Os03g03100 (OsMADS50, Supplementary Figure 4). It was noted that haplotype 2 (pink group) had the highest number of accessions, with 66.5%, followed by haplotype 1 (brown group) with 18.0% and finally by haplotype 3 (blue group) with 15.5%. Nipponbare and Sherpa clustered with haplotype 3 along with 60.4% of temperate japonica (tempjap) present in the 3KRGP, while haplotype 1 had 39.2% of the tempjaps (as well as IRAT109) and haplotype 2 only had one tempjap accession.

The current paper evaluated RILs derived from the LT tolerant Australian rice variety, Sherpa for variation in SS and other key traits in controlled temperature glasshouse conditions. Evaluation of the full set of RILs resulted in a mean SS of 54.7 and 37.9% from Experiments 2 and 4, respectively and similarly, using a set of differential genotypes revealed a mean SS of 51.0% (Experiment 6). Across the two full population evaluations, a total of 13 RILs were identified to perform well with mean SS below 30%. The mean SS from our current experiments were comparable with previous results (Susanti et al., 2019) using the similar screening method developed by Mitchell et al. (2016). Likewise, the heritability values of SS from the experiments conducted in this study ranged from 0.72 to 0.90 and were also similar to a previous report (Susanti et al., 2019). This consistently high heritability of the trait along with the existence of large genetic variation (Guo et al., 2020) indicates that SS may have a good response to selection. Additionally, Ha et al. (2018) were able to show consistency in low-temperature tolerance between AP and PW indicating that current results can be utilized in both growing conditions. To facilitate the efficient development of these genotypes, it is imperative to identify genomic regions controlling SS for their use in molecular marker-assisted breeding programs.

To our knowledge, for the first time, the current study has shown the utility of an impedance flow cytometer in measuring TP of rice both in the context of ideal and low-temperature conditions. Additionally, tolerant genotypes with low SS showed higher NDA and TP, indicating that these two floral traits may be used as surrogates to select for genotypes with low SS in low-temperature conditions. NDA is relatively easy to phenotype and would require a simple stereomicroscope while TP measured using a flow cytometer (Heidmann et al., 2016) can be collected in minutes, both several weeks before the data for SS after panicle maturity becomes available. We utilized a set of differential genotypes to dissect the relationships between floral traits and the development of low SS. In the experiment where genotypes were exposed to LT, it was shown that the low SS tolerant group developed a higher number of pollen per spikelet resulting in higher NDA and eventually higher POS. The importance of pollen number in anther in determining SS after exposure to LT conditions has been demonstrated in previous studies (Mitchell et al., 2016; Susanti et al., 2019) and was supported by the fact that pollen swelling was found to be the prime force for anther dehiscence (Matsui et al., 1999). It is important to note though that this relationship with SS was absent in warm ideal growing conditions indicating that the previous relationships observed were only induced after exposure to LT conditions. Farrell et al. (2006b) also noted that SS was related to a number of engorged pollen and another length only after cold water exposure. Significant correlation for TP between warm and LT conditions also existed although there was no significant association between TP from warm and SS of LT experiment. This suggested that the genotypic variation in TP in warm conditions affected the variation in TP in LT treatment to some extent but was not strong enough to affect SS in LT conditions and ultimately, LT tolerance.

Three genomic regions controlling LT tolerance at YMS were identified, located at chromosomes 3 (qYMCT3), 8 (qYMCT8), and 4 (qYMCT4), and were found to co-locate with previously reported QTL. Several genomic regions associated with LT tolerance at various stages of crop development have been previously identified. In terms of the reproductive stage, a few regions have been identified to be associated with tolerance at either YMS or the booting stage (Ye et al., 2010; Lei et al., 2018; Proud et al., 2022). QTL qYMCT3 was found to be located in the short arm of chromosome 3 roughly translating to between 1.25 and 2.98 Mb physical map position. This QTL was found to co-locate with qCTSF3.1, a genomic region associated with LT tolerance at the booting stage identified using chromosome segment substitution lines derived from a cross between japonica and indica varieties (Lei et al., 2018). By utilizing a genome-wide association mapping approach, Guo et al. (2020) identified SNPs located in 15.87–19.49 Mb of chromosome 8 using 522 diverse cultivars. Finally, the confidence interval of qYMCT4 was also found to co-locate with a genomic region identified by examining 168 diverse accessions in high altitudes (Guo et al., 2020). qYMCT4 was identified only in Experiment 2 and may suggest that this region is significant in cooler/harsher temperatures. A similar observation was made in DTH QTL identified which may signify the effect of these genomic regions in the short days (Experiment 3), in long day/natural light conditions (Experiment 1), or both. Identification of these inherent useful genomic regions from Sherpa should be helpful in (1) identifying existing varieties with favorable allele or QTL combinations and (2) introgressing additional LT tolerance genes/QTL to further improve the tolerance levels of the variety or subsequent development of improved genotypes with multiple LT tolerance QTL. Pyramiding of favorable alleles of high-value genomic regions and cloned genes (CTB4a, Ctb1, bZIP73, and OsAPX1) have been shown to increase LT tolerance of genotypes in North China (Guo et al., 2020). Transgressive segregation in both DTH and SS existed. Specifically, in terms of SS, both parents contributed genomic regions and alleles associated with low SS, which may explain the transgressive segregation observed in this population. Although rare, transgressive segregation for SS after LT exposure at YMS or booting stage has been previously observed (Takeuchi et al., 2001). The final QTL model for SS in the current study explained around a quarter of the phenotypic variation in the population. The construction of specialized mapping populations such as chromosome segment substitution lines may aid in the identification of additional genomic regions governing SS (Yamamoto et al., 2009).

There is evidence of a pleiotropic effect of the QTL/gene controlling SS for LT tolerance (qYMCT3) and DTH (qDTHW3.1) in this population. To our knowledge, our subsequent analyses are the first attempt in linking OsMADS50 with LT tolerance in rice. Several efforts using genomics, transcriptomics, and gene functional analyses conducted in rice and Arabidopsis thaliana were in congruence with the observations made in our current study. Using a set of genotypes derived from the IRRI breeding program, qDTHW3.1 was also found to be associated with DTH in the indica subpopulation (Begum et al., 2015). The same region was also identified using a NIL derived from Oryza glaberrima (Bian et al., 2011). Additionally, it was pointed out that LOC_Os03g03100 was the best candidate for this region (Begum et al., 2015). WGRS comparative analysis in our current study supported this claim by showing the absence of polymorphism between parents in Ehd4, while several polymorphisms were identified in LOC_Os03g03100. The primary candidate gene for qYMCT3/qDTHW3.1 was therefore LOC_Os03g03100 encoding for a MADS-box transcription factor (OsMADS50). Functional analysis of OsMADS50 showed that it is an ortholog of A. thaliana’s SUPPRESSOR OF OVEREXPRESSION OF CO1/AGAMOUS-LIKE 20 (SOC1/AGL20) (Lee et al., 2004). Their experiments also showed that overexpression of OsMADS50 resulted in early flowering while RNAi and knock-out plants exhibited phenotypes with late flowering (Lee et al., 2004; Ryu et al., 2009). MADS-box gene family has been identified as important molecular players in abiotic stress responses such as LT (Castelan-Munoz et al., 2019). SOC1 is a known key floral activator integrating multiple floral inductive pathways in A. thaliana and was in parallel to the observation in rice through its ortholog, OsMADS50. SOC1 was also shown to participate in LT sensing. Loss of function analysis (soc1-2) showed induced expression of LT responsive genes such as cold regulated and C-repeat/dehydration response element binding factors, while its overexpression resulted in decreased expression of such genes, indicating that SOC1 was a negative regulator of the expression of LT response genes (Seo et al., 2009). The above series of experiments demonstrated the crosstalk between flowering and LT sensing in A. thaliana, and with the shared genomic regions identified in the current study for DTH and SS, it may very well be the same case in rice, for OsMADS50. Co-expression and GO enrichment analyses of genes associated with the expression of OsMADS50 showed its link with calcium ion transport. Cold signal transduction has also been shown to involve the activation of calcium ion channels and pumps; and increasing evidence shows that calcium signaling plays a crucial role in conferring LT tolerance among plants (Yuan et al., 2018), further strengthening the plausible crosstalk between flowering and LT tolerance as mediated by OsMADS50. However, further experimentation utilizing over-expression and knock-out mutants of OsMADS50 evaluated in LT environments is needed to underpin this crosstalk and relationship in rice. Additionally, while the current analysis and literature point out OsMADS50 as the most prudent candidate gene for qYMCT3 and qDTHW3.1, it is still also possible that another gene or independent genes contribute to LT tolerance and DTH in this population.

With these genomic regions identified, their validation in multiple genetic backgrounds is an imperative step for their utility in respective breeding programs (Cobb et al., 2019). As the plausible useful allele OsMADS50 may result in delayed heading, it is also important to identify additional genomic regions which may nullify this effect. Although late to head, osmads50 mutant plants did eventually flower indicating that other positive regulators induced flowering in the absence of the expression of OsMADS50 (Ryu et al., 2009). This is in congruence with our data which identified other genomic regions associated with earlier flowering from Sherpa. The presence of the other five genomic regions with Sherpa having the favorable (early flowering) allele suggested that other genomic regions may be used to produce LT tolerant, early flowering genotypes from this donor. Furthermore, haplotype analysis showed the possible utility of OsMADS50 favorable haplotype across rice subpopulations. The Sherpa haplotype (haplotype 3) was shown to have a relatively high frequency (60%) in the tempjap subpopulation suggesting that although it might be useful to introgress this haplotype into genotypes with non-favorable haplotype, the utility of this gene to improve LT tolerance in tempjaps will be of higher magnitude in the survey and identification of genotypes already possessing this useful haplotype. This type of utility was highlighted recently in an elite breeding core panel where several high-value disease and insect resistance genes were present at considerable frequencies (although not fixed) and thus represent haplotypes/alleles that are easily selectable in existing breeding populations (Juma et al., 2021). After successful QTL validation, molecular markers could be designed for their survey and introgression in breeding programs (Vinarao et al., 2021a). An assessment of high throughput molecular markers in rice showed that only two QTL for seedling LT tolerance, SCT1, and COLD1, were developed for this purpose³ indicating the opportunity to develop markers for reproductive stage LT tolerance. Therefore, the development of high throughput molecular markers for the three QTL identified in this study as well as known cloned LT tolerance genes should aid in the tailored genotype identification and pyramiding of additional beneficial genomic regions to develop AP adapted, LT tolerant genotypes.