Breeding for Low Temperature Germinability in Temperate Japonica Rice Varieties: Analysis of Candidate Genes in Associated QTLs

: In temperate areas, rice deals with low temperatures that can affect plant growth and crop yield. Rapid germination is required for adequate plant establishment in the field, therefore obtain-ing cultivars that maintain this phenotype under suboptimal temperature conditions is a challenge for rice breeders. Our study aimed to investigate temperature-induced expression changes in genes underlying quantitative trait loci (QTLs) associated to this trait (low temperature germinability, LTG) that were detected in a previous genome wide association study (GWAS). In the context of a breeding program for japonica rice cultivars adapted to cultivation in Spain, we obtained two biparental families of lines derived from hybridization with two cold tolerant Italian cultivars, and we have studied the effect on the LTG phenotype of introgressing these QTLs. A wide region in chromosome 3 was related to significant increases in seedling growth rate at 15 °C, although the extent of the effect depended on the analyzed family. In parallel, we studied the pattern of expression during germination at different temperatures of 10 genes located in the LTG-associated QTLs, in five japonica rice cultivars and in a biparental family of recombinant inbred lines (RILs). Cold induced changes in the expression of the 10 analyzed genes, with significant differences among genotypes. Variation in LTG phenotype was consistently associated with changes in the pattern of expression of five genes from the tagged regions in rice chromosome 3, which encoded for enzymes implicated in phytohormone metabolism (OsFBK12, Os3Bglu6), oxidative stress (SPL35, OsSRO1c) and Mn homeostasis maintenance (OsMTP8.1). Differential expression induced by cold in two regulatory genes (Os02g0824000 and Os06g06400) also contributed to explain low temperature tolerance during rice germination. In conclusion, introgression in defective cultivars of favorable alleles for these genes would contribute to the genetic improvement of LTG in japonica rice varieties.


Introduction
Different cultivation systems and plant adaptation strategies have allowed rice (Oryza sativa L.) crops to expand worldwide, from its tropical region of origin in Southern Asia to temperate zones in North America and Europe. Nowadays, 90% of global production comes from Asia and 75% is harvested in irrigated fields [1], this cultivation system is highly demanding on water resources. In the European Union, rice is cultivated from about 450,000 ha, mainly in Italy and Spain, but also in Greece, Portugal, and France, with a mean yield of 7 t/ha. In these countries, rice crops are usually located in soils with salinity and flooding problems, in or near wetlands under high environmental protection.
About 80% of the rice cultivated area in the EU are cultivars belonging to the japonica group of O. sativa. Japonica rice cultivars differ from indica types in terms of plant architecture, specific agroecological adaptations, and grain quality traits [2], and have expanded from tropical to temperate areas during domestication and breeding processes, despite the narrow genetic pool [3]. Temperate japonica varieties progressively adapted to colder climates with shorter growing seasons by flowering under long day conditions [4].
Breeding rice for EU producers is an activity conditioned by several factors, such as consumers' demand of specific grain types, the changes in the Community Agricultural Policy (CAP), the adoption of integrated management practices, the demand for reduction of water consumption of the crop, and, more recently, the introduction of both the Clearfield ® technology to control red rice and hybrid cultivars. Therefore, the specific objectives of the developed breeding programs are to improve yield, plant adaptation, resistance to pest, and diseases, and grain quality must meet these requirements.
Taking this into consideration, in Europe, as in other temperate areas worldwide, low temperature is a major limiting factor for rice productivity [5]. This plant species is more sensitive to cold stress than other cereals, particularly when temperature drops below 18-20 °C in critical phases, such as during germination, crop establishment, and the reproductive stage. Cold stress affects seed germination and seedling establishment, which may reduce the growth of young plants and crop yield to some extent, since good establishment is a key factor for the success of subsequent crop growth [6]. The genetic improvement of low temperature tolerance at germination (LTG) is therefore an objective of rice breeding programs in these regions. Cold tolerance is a heritable trait for which genetic variation has been reported to be associated with the geographical distribution of rice cultivars, being that those from the japonica group are more tolerant than indicas [7]. Studies on rice cold tolerance have employed a variety of methods to evaluate this trait, mainly depending on the stage of development of the plant [6]. At the germination stage, the germination vigor and the seedling survival rate are the two main criteria used to evaluate cold tolerance [8]. According to these criteria, some European cultivars, such as Italica Livorno [9] and Maratelli [10], have been described as cold tolerant and used as donor parent to introgress the trait in defective plant material. More recently, an accession of wild rice (Oryza rufipogon) has been also proposed as potentially useful donor for LTG improvement [11].
When breeding rice for traits, such as seed survival rate and plumule growth or greening under low temperature conditions, moderate to high narrow and broad sense heritability estimates have been reported and several genes and QTLs (quantitative trait loci) have also been associated with cold tolerance of rice seedlings [1,6]. Although rice seed germination ability mainly depends on cultivar heredity, both exogenous and endogenous factors, including water availability, temperature, light, circadian rhythm, and phytohormones, participate in its regulation [12]. Most of these mapping studies, which have detected QTLs regulating this trait across the 12 rice chromosomes, were performed using biparental populations derived from indica × japonica crosses, but recently genome wide association studies (GWAS) using different collections of rice cultivars have also identified many chromosome regions regulating this complex character [13][14][15][16][17][18][19][20]. However, candidate gene identification in these regions is hampered by a limited mapping resolution, since rice is an autogamous plant, and therefore has a modest rate of decay of linkage disequilibrium (~200 kb), which usually contains 10-20 genes [21]. For example, two glutathione S-transferases were reported as the candidate genes underlying qCTS12, a QTL regulating cold tolerance of rice seedlings [22]. Expression analyses allowed the identification of candidate genes in three QTLs associated to LTG: a protein with unknown function in qLTG-3-1 from Italica Livorno probably involved in tissue weakening during germination, which is expressed in the initial phases of the process [23]; a zinc finger protein in qLVG7-2 [18], and both a glycoside hydrolase and an NBS-LRR domain containing protein in qLTG_sRDP2-10a [24]. Cold-stress related genes have been cloned from the rice genome, including COLD1, which encodes a regulator of G-protein which activates the Ca 2+ channel for sensing low temperature and to accelerate G-protein GTPase activity [25], and the gene OsSRFP1, that may have dual functions in modulating abiotic stress responses in rice, at least in part, by negatively regulating antioxidant enzymes that remove reactive oxygen species [26]. Understanding the complex molecular regulation of the response of rice seeds to low temperature at germination is therefore a challenge for geneticists and breeders. The current knowledge about QTLs regulating rice cold tolerance and the implementation of genomic breeding and gene/QTL pyramiding strategies for the genetic improvement of this trait has been recently reviewed [27].
In our previous work [16], we used a collection of 180 temperate japonica rice varieties and a panel of 1672 SNP markers to perform a GWAS that identified 24 loci associated to higher growth rates in seedlings germinating at 15 °C. Among them, a core panel of 13 SNP markers consistently associated with rice LTG were detected in chromosomes 2, 3 (5 markers), 6 (4 markers), 7, 10, and 11. The haplotype of the cultivar Italica Livorno for this set of SNPs was partially shared by other Italian cold tolerant cultivars, such as Carnise, Maratelli, Nuovo Maratelli, and Poseidone, and the number of favorable alleles correlated with LTG. We proposed these markers for developing more adequate breeding programs, particularly for japonica rice cultivars, since introgression by backcrossing qLTG-3-1 from Italica Livorno or Maratelli has not been effective, and therefore pyramiding several QTLs is probably needed in order to significantly increase LTG [10,28].
Aiming to better characterize the detected positive alleles, we also examined the 200 kb region centered on the 24 tagged SNPs looking for potential causal genes, and found some encoding enzymes involved in mechanisms of stress tolerance (thioredoxins and leucine rich repeat proteins), and in controlling development (pentatricopeptide repeat domain containing proteins or transcription factors), as well as transposable elements. Furthermore, three of the tagged gene products (Os02g57800, Os03g12820, and Os06g06400) had been previously reported as being differentially expressed in stressed rice seedlings [29].
In the present work, we have selected ten of these genes for performing expression studies using japonica rice cultivars with contrasting LTG responses: three tolerant Italian cultivars, Italica Livorno, Carnise and Nuovo Maratelli, the cultivar Guadiamar, a coldsensitive variety cultivated in Northeast Spain, as well as the Spanish variety Montsianell, which is commonly used by growers to avoid germination problems because of its good establishment in the field. This cultivar has no Italica Livorno haplotype for the marker core panel, and only presents an additional favorable allele (that in chromosome 7) than the cultivar Guadiamar. To validate results from our GWAS, we also analyzed an advanced recombinant family (F6) derived from the cross between Guadiamar and Italica Livorno. These cultivars differed by 10 SNPs from the core panel, those in chromosomes 3, 6 and 7. The core panel was also employed to characterize F3 families derived from the hybridization between other two cultivars, Copsemar 7 and Nuovo Maratelli, which differed in 7 markers.

Selected Rice Genome Regions and Candidate Genes
From the GWAS results, performed to identify rice genome regions associated with LTG [16], we selected five temperate japonica cultivars to analyze the expression of the detected candidate genes: four cultivars that had shown some degree of tolerance (Carnise, Italica Livorno, Montsianell and Poseidone), and one cold sensitive variety (Guadiamar) with slow germination. Table 1 shows the haplotype of these cultivars for 15 SNPs markers significantly associated to LTG in japonica rice, being 13 of them the core panel with stronger association. The table also shows the ten investigated candidate genes in their tagged region. We tried to study other genes located in tagged regions from chromosomes 6 and 7, but amplifications failed to give specific products as was expected.

Gene Expression Analyses in Rice Varieties
Germination experiments were carried out, as described in [16], using rice seeds provided by COPSEMAR (Sueca, Valencia, Spain). For each cultivar, we determined mean seedlings length at 15 °C and at 25 °C, after 14 and 21 d, and after 4 and 6 d, respectively. We estimated growth rates at different temperatures (named V15 and V25) as the ratio seedling length/number of days, using three replicates with 25 seeds each. For comparative purposes, other plates were also tested at optimal temperature (30 °C), being monitored for 4 days.
For each rice cultivar, RNA samples were obtained using the Plant and Fungus RNA isolation kit (Norgen Technologies, Thorold, ON, Canada) from 30 mg of seedlings without roots, growing either at 25 °C (4 and 6 days after sowing), and at 15 °C (14 and 21 days after sowing). For the cultivars Guadiamar and Poseidone it was possible to collect enough material only at the second sampling time. RNA samples were treated with DNase (Turbo DNA-free kit, Ambion, Life Technologies, Thermo Fisher Scientific, Waltham, MA, USA) and then cDNA synthesis was performed using Transcriptor First Strand cDNA Synthesis kit (Roche Life Science, Germany). Both processes were performed following manufacturer's specifications.
For real time PCR amplifications, primers were designed using the sequences found in the Rice Genome Annotation Project database (http://rice.plantbiology.msu.edu/, accessed February 16, 2017) for the candidate genes (Supplementary Table S1). Two rice genes were tested as internal reference, the transcription factor 1α and the 18S RNA, using primers previously described [30]. The efficiency and specificity of the all primer pairs were evaluated. Amplifications were performed in a Step One Plus Real time PCR system (Life Technologies, Thermo Fisher Scientific, Waltham, MA, USA) using the SYBR Pre Mix Ex Taq II kit (Takara Bio, Kusatsu, Shiga, Japan) according to the manufacturer's recommended protocol with three technical replicates, and using a program consisting of an initial denaturalization step at 95 °C × 30 s, following by 40 cycles of 95 °C × 5 s and 60 °C for 30 s. Primer efficiencies were always greater than 88%. Gene expression levels were estimated by the 2 −∆CT method [31] relative to threshold cycle (CT) values of the gene coding 18S RNA. Differential expression rates were estimated as the base 2 logarithm of the ratio between 2 −∆CT values at 15 °C and at 25 °C (that is, Log2 of fold change).

Phenotypic and Genetic Characterization of the Guadiamar × Italica Livorno Family
Petri dishes with 25 seeds each of 73 F6 Recombinant Inbred Lines (RILs) derived from the cross Guadiamar × Italica Livorno, were prepared and incubated at 15 °C or 30 °C (3 replicates each). Seedlings' lengths were recorded after 14 and 21 days at 15 °C, and after 4 days at 30 °C, respectively, to estimate growth rates as seedling length/number of days (V15 and V30, respectively). Averaged values were obtained for each line for the complete assays and for the third week of the assay at 15 °C (V153w). We also estimated mean percentage of growth rate at 15 °C as compared to that at 30 °C (%V15/V30).
For genotyping these RILs, DNA from each line and from the two parental lines was isolated using leaves sampled in 2-week old plants grown in peat-moss:vermiculite (1:1) substrate at 30 °C. A modified CTAB protocol routinely used in our lab was deployed: about 100 mg of leaf tissue was frozen with liquid nitrogen and then ground in a mortar using a pestle; the obtained powder was transferred to a 2 mL tube containing 700 μL of extraction buffer (2% CTAB, 2% PVP, 0.1 M Tris-HCl, 0.2 mM Na-EDTA and 1.4 M NaCl), 50 μL of lysis buffer (5% N-lauryl sarcosine, 0.1M Tris-HCl) and 5 μL of a 20 mg/mL Proteinase K solution (Norgen Technologies, Canada). After vortexed for 30 s, samples were incubated at 60 °C for 30 min and then centrifuged (12,000 rpm × 10 min). The supernatant was recovered, transferred to a new tube, and then extracted with 470 μL of chloroform:isoamyl-alcohol (24:1). After a new centrifugation step, the aqueous phase was transferred to a new tube, and nucleic acids precipitated by adding 0.7 V of isopropanol and recovered by centrifugation (12,000 rpm × 10 min, 4 °C). Each sample was diluted in 50 μL of ultrapure water and treated with RNase. The Sequenom MassARRAY technology (SCAIE, University of Valencia, Spain) was employed to determine the genotype of each F6 line for a panel of 10 SNPs that differed between both parents (Guadiamar and Italica Livorno), from the core panel of 13 markers significantly associated with LTG and also for a SNP at the position 28,375,652 in the chromosome 7, because a gene in this region was included in the expression analyses (Table 1).
Finally, 16 RNA samples were prepared as described previously, using 21-day-old rice seedlings grown at 15 °C, and 4-day-old seedlings grown at 30 °C, from the two parents and from 6 RILs (6, 11, 18, 21, 37 and 40) selected according to their genotype for the core panel. cDNA synthesis and real-time amplifications were performed following the same protocols used in the experiments with 5 rice cultivars.

Selection for LTG in Copsemar 7 × Nuovo Maratelli F3 Families
A further germination experiment was performed by incubating at 15 °C Petri dishes containing 25 seeds of each 13 F2 plants derived from the cross Copsemar 7 × Nuovo Maratelli (325 seeds in total). After 21 days, germinated seedlings were transferred to peatmoss:vermiculite (1:1) substrate and maintained at 30 °C for two weeks before isolating DNA from leaves, as described before. The Sequenom MassARRAY technology was therefore used to genotype 86 F3 plants for the 7 SNP markers, located in chromosomes 3 and 6, which showed different alleles in the parental cultivars (Table 1). Plants were finally transferred to the greenhouse and cultivated to obtain seed of F4 families, which will be further used in our breeding scheme.

Statistical Analyses
Seedling growth rates, as well as levels of gene expression (2 −ΔCT ) and Log2-fold change values in experiments performed at different temperatures, were analyzed by ANOVA (normal data) or by non-parametrical tests (non-normal data) to determine the existence of significant differences (α = 0.05), using IBM Statistics SPSS software. We also analyzed correlations (Spearman coefficient, ρ) among seedling growth parameters.

Two Temperatures Phenotype Analysis in Rice Varieties
Seedlings of the five tested rice varieties differed in growth rate depending greatly on temperature, as expected from our previous studies [16]. Seedling growth rates decreased when germination was carried out at suboptimal temperatures, 15° or 25 °C, particularly at 15 °C, as compared to those observed in rice seeds germinating at 30 °C (Supplementary Figure S1). We observed significant differences among cultivars (p < 0.001) for seedling length after 4 and 6 days at 25 °C, and after 14 and 21 days in the experiment performed at 15 °C ( Figure 1). Best results at 15 °C were obtained for Carnise and Italica Livorno cultivars after 14 days (5 ± 1 mm) and 21 days (11 ± 2 mm), while Montsianell and Poseidone showed intermediate response (4 ± 1 and 8 ± 2 mm after 14 and 21 days, respectively). In the experiment performed at 25 °C, significantly higher values for seedling length were determined in seeds from Carnise, Italica Livorno, and Montsianell; therefore, germinating vigor was superior in these three varieties. The cultivar Guadiamar showed the poorest results with the lowest values of seedling length at both temperatures ( Figure  1). However, Poseidone was competitive with Montsianell in seedling length at 15 °C.

Gene Expression Analyses in Rice Varieties
After real time amplifications of cDNA samples prepared from seedlings of the five rice cultivars germinating at either 15 °C or 25 °C, we observed that the expression of the 1α elongation factor gene was affected by temperature, and therefore only 18S RNA gene results were used as internal reference to estimate 2 −∆CT values. The ten studied genes were expressed during germination at both temperatures, with higher expression levels being observed on average for Os02g57800, Os03g07530, Os03g11420, Os03g12530, Os03g12820, and Os11g05540. Gene expression was significantly affected by temperature, as well as by the sampling time, in most of the cases ( Table 2).   (Table 2), seedlings of Montsianell, Carnise, and Italica Livorno showed significantly higher expression levels when germinating at 15 °C than at 25 °C for two genes located in chromosome 3, Os03g12530 and Os03g12820 (Figure 2a), and for Os11g05540 (Figure 2b), while the remaining genes showed significant increases for only one or two of these cultivars. Differential expression rates varied among cultivars and genes, being on average higher for Carnise (Log2-fold change values were 1.34 ± 0.03, 2.12 ± 0.10, and 1.51 ± 0.00 for Os03g12530, Os03g12820, and Os11g05540, respectively) than for Italica Livorno (0.43 ± 0.12, 2.09 ± 0.03, and 0.46 ± 0.21), and for Montsianell (0.65 ± 0.37, 1.72 ± 0.18, and 0.34 ± 0.18). At the second sampling time, after 21 days at 15 °C and after 6 days at 25 °C, in the five studied cultivars ( Table 2) we observed higher expression levels at 15 °C than at 25 °C for three genes: Os03g12820, Os06g06400 and Os07g47470 (Figure 3). However, differential expression rates varied among cultivars, since for the Os03g12820 gene we estimated for the cultivar Guadiamar a Log2-fold change of 1.93 ± 0.09 while the Italian cultivar Poseidone reached a value of 3.55 ± 0.15. Differential expression rates of the gene Os06g06400 ranged from 0.65 ± 0.29 in Carnise to 4.30 ± 0.46 in Italica Livorno, while Log2fold change estimates for Os07g47470 varied from 1.20 ± 0.26 in Guadiamar to 4.78 ± 0.12 in Italica Livorno. Interestingly, five genes showed opposite expression patterns between the four cultivars exhibiting LTG and the cultivar Guadiamar (Figure 3). The three Italian cultivars and Montsianell showed significantly higher expression levels under cold conditions than at 25 °C for the other four genes located in chromosome 3, and for the analyzed gene in chromosome 11, being the highest Log2-fold change values observed for Italica Livorno seedlings: 4.27 ± 0.18 in Os03g07530, 3.73 ± 0.44 in Os03g10750, 3.94 ± 0.12 in Os03g11420, 3.75 ± 0.15 in Os03g12530, and 5.28 ± 0.09 in Os11g05540. In Guadiamar seedlings, when germination was carried out at 15 °C these genes either remained unchanged, as observed for Os03g10750, Os03g11420, and Os03g11420, or were significantly repressed: Log2-fold change = −0.90 ± 0.20 for Os03g7530, and Log2-fold change = −1.02 ± 0.13 for Os03g12530. Finally, at this sampling time, genes Os02g57800 and Os10g40390 were overexpressed under cold conditions in Italica Livorno, Montsianell and Poseidone, but not in Carnise (Figure 3).

LTG in the Guadiamar × Italica Livorno F6 Family
Seedlings from the 73 RILs of the Guadiamar × Italica Livorno family differed significantly (p < 0.001) in their growth rate in germination experiments performed at either 30 °C (V30; optimal conditions) or 15 °C (V15; cold conditions). We estimated V30 values that ranged from 1.51 ± 0.12 to 3.74 ± 0.22 mm/d, with an average of 2.52 ± 0.34 mm/d ( Figure  4a), while V15 values varied from 0.16 ± 0.02 to 0.59 ± 0.07 mm/d, with an average of 0.33 ± 10 mm/d (Figure 4b). Growth rates at 15 °C positively correlated with growth rates at 30 °C (ρ = 0.336, p = 0.004), higher V15 rates being explained by faster seedling development during the third week of the experiment, since we also estimated a significant positive correlation with V153w values (ρ = 0.854, p < 0.001). Significant differences were also observed among lines for these values (p < 0.001), and for the percentage of seedling growth rate maintained under cold conditions (p < 0.001), since this estimate (%V15/V30) correlated directly with V15 data (ρ = 0.629, p < 0.001) and inversely correlated to V30 data (ρ = −0.455, p = 0.004).  Significant effects of the genotype on mean seedling growth rates of the 73 RILs were only detected for two out of 12 SNP markers. Seedling growth rates of lines with the Italica Livorno allele for the SNP located at position 3,828,648 in chromosome 3 were less affected by cold than those with the Guadiamar allele (p = 0.005), since we estimated them to have higher %V15/V30 values (Figure 5a). This result was explained by higher V15 values observed in RILs carrying the Italica Livorno allele, although when compared to lines with the Guadiamar allele for this SNP, differences were not significant for this parameter (p = 0.062). In contrast, lines with the Italica Livorno allele at position 5,506,710 in chromosome 3 showed significantly higher values of V15 (Figure 5b) than those that inherited the Guadiamar allele (mean 0.33 ± 0.10 front to 0.27 ± 0.04 mm/d, p = 0.027), as well as higher V153w rates (p = 0.004), and %V15/V30 (p = 0.008) estimates (data not shown). Finally, gene expression analyses were performed using cDNA samples prepared from seedlings germinating at either 15° or 30 °C after 21 and 4 days, respectively. We analyzed differential expression of the selected 10 genes in both parents, and in six F6 lines: three lines with the Guadiamar haplotype for the panel of 10 SNPs markers (6 and 11), or for 9 out of these markers (line 18 carries the Italica Livorno allele for the SNP at position 13,128,278 in chromosome 7), and three lines (21, 37 and 40) with the Italica Livorno allele for the 10 SNPs. These two groups showed significantly different values of seedling growth rates at both temperatures of germination: mean V15 = 0.26 ± 0.04 front to 0.52 ± 0.11 mm/d (Supplementary Figure S2), and mean V30 = 1.78 ± 0.21 front to 2.93 ± 0.44 mm/d, respectively. Although Log2-fold change rates varied among lines, we found significant differences between both groups of RILs with contrasting LTG phenotypes for the ten analyzed genes ( Figure 6). When germination was carried out at 15 °C, genes Os03g07530, Os03g10750, Os03g11420, and Os06g06400 were significantly overexpressed in F6 lines with the Italica Livorno haplotype, while these four genes were significantly repressed in lines with Guadiamar alleles. Cold also significantly increased gene expression of Os03g12530 and Os03g12820 in seedlings of the six analyzed lines, but fold change was significantly greater in the group with the Italica Livorno alleles than in that presenting the Guadiamar alleles (on average Log2-fold change was 2.20 ± 0.56 front to 0.61 ± 0.98, p < 0.001, and 3.09 ± 0.76 front to 1.33 ± 0.91, p < 0.001, respectively). In contrast, under cold conditions expression of the Os02g57800 gene was significantly reduced when compared to germination at optimal temperature, and this decrease was greater in lines carrying the Guadiamar haplotype than in those with the Italica Livorno haplotype (Log2fold change mean values of −2.57 ± 1.22 front to −0.57 ± 0.58, p = 0.001). Cold also reduced expression of genes Os07g47470, Os10g40390, and Os11g05540 in F6 lines with the Italica Livorno haplotype (Figure 6), while the opposite trend was observed for lines with the Guadiamar haplotype, in which expression of these genes remained on average unaffected (Os10g40390 and Os11g05540) or was significantly increased (Os07g47470). Figure 6. Differential expression (Log2 of fold change) of 10 rice genes in 21-day-old seedlings germinating at 15 °C, and in 4-day-old seedlings germinating at 30 °C. Analyses were performed using seedlings of the japonica rice cultivars Guadiamar (GU) and Italica Livorno (IL), and from 6 F6 lines derived from its cross. Lines 6, 11 and 18 shared the GU haplotype for a panel of 10 SNP markers associated with LTG, and lines 21, 37, and 40 presented the IL haplotype. Data are expressed as mean ± SD of gene expression estimates derived from amplifications performed in triplicate. * indicates the existence of a significant difference between both groups of lines.

Selection for LTG in Copsemar 7 × Nuovo Maratelli F3 Families
Analyses of LTG were also performed with 325 seeds from 13 plants selected among the F2 family derived from the cross Copsemar 7 × Nuovo Maratelli. Parental lines differed in seedling growth rate under cold conditions, that was on average 0. 36 (5,874,917 and 6,639,548) in chromosome 3 showed significantly higher seedling growth rates at 15 °C (V15) than 11 F3 plants that inherited at both loci the Copsemar 7 allele (p = 0.027 and p = 0.013, respectively). Most of these plants with the LTG phenotype also presented the Nuovo Maratelli allele at position 5,506,710 (52 out of 57). Although results might be biased by the phenotypic selection performed, we found that plants with Nuovo Maratelli alleles for these 3 SNPs showed higher V15 rates than 10 plants that inherited the Copsemar 7 haplotype (Figure 7), since seedling growth rates of plants increased on average from 0.27 ± 0.12 to 0.40 ± 0.16 mm/d.

Discussion
Natural variation in cold tolerance of rice cultivars results from plant adaptation out from its tropical region of origin. Cold acclimation response depends mainly on biosynthesis of osmotics, generation of antioxidants, increase of membrane fluidity, and maintenance of cell wall integrity. The physiological changes in rice under cold stress were reviewed in [19], and this abiotic factor usually promoted the greatest changes in the overall metabolism as compared to salinity and iron toxicity [32]. Low temperature signal is transduced by phosphorylation and transcriptional cascades, which also implicate ethylene biosynthesis, resulting in a large re-programming of the plant transcriptome [33]. Moreover, cold acclimation is a very complex response that is actually regulated at all levels, i.e., transcriptionally, post-transcriptionally, translationally, and post-translationally. Therefore, the trait is subjected to both genetic and epigenetic control. Regulation of gene expression and signaling transduction under cold stress is mediated by transcription factors from several families, and also by proteins involved in RNA processing and nuclear export. Two different regulatory pathways have been proposed to complement each other and act sequentially to adapt rice to immediate and persistent cold stress: response of plants to cold stress is mainly activated by CBF-DREB transcription factors (C-repeat binding factor-dehydration responsive element binding) that are most likely crucial in responding to short-term cold stress, and a MYBS3 transcription factor-mediated system regulated by abscisic acid (ABA) that is more important for long-term adaptation to persistent cold stress [34]. Furthermore, studies performed with japonica cultivars suggested that during cold conditions seedling growth is inhibited by increased ABA levels, and therefore cold tolerance was associated with maintaining relatively low levels of ABA [35].
Interestingly, other authors [36] reported that different mechanisms of gene regulation can be operating in the two subspecies of rice. Thus, the tolerance to cold stress in indica rice was found to be regulated via the CBF-dependent as well as the CBF-independent pathway, while cold tolerance in japonica occurs primarily via the CBF-independent pathway, being mediated by ABA.
In our experiments, we first used seeds of five japonica rice cultivars with different LTG phenotypes to study patterns of gene expression in seedlings germinating under cold conditions (15 °C), as compared to control conditions (25 °C). Ten genes were selected in regions tagged as associated to LTG in a previous GWAS [16], which included some Italian cultivars that exhibit LTG (Carnise, Italica Livorno, Maratelli, Nuovo Maratelli, and Poseidone). On average, these 10 genes were upregulated under cold conditions at least during seedling elongation (after 21 days at 15 °C as compared to 6 days at 25 °C). Quantitative RT-PCR analyses allowed the identification of candidate genes underlying mapped QTLs controlling cold tolerance-related traits in rice [24], and to dissect mechanisms for cold tolerance operating in japonica and indica subspecies [36]. Reduction in seedling length after cold stress and differential expression analyses also discriminated cold tolerant cultivars among the USDA rice mini-core collection [37]. In our study, higher expression levels at 15 °C were detected in four genes located in chromosome 3: Os03g07530, Os03g11420, Os03g12530 and Os03g12820, and in the gene Os11g05540 at chromosome 11, and, particularly, in seedlings of the reference cold tolerant cultivar, Italica Livorno.
Four of these genes, Os03g07530, Os03g11420, Os03g12530, and Os11g05540, as well as Os03g10750, showed opposite trends of differential expression between four cold tolerant cultivars (Carnise, Italica Livorno, Montsianell and Poseidone), in which were upregulated under cold conditions, and the cultivar Guadiamar, for which expression levels remained unchanged or were repressed (Figure 3). This pattern of expression was partially corroborated in the experiments performed with Guadiamar × Italica Livorno RILs, since under cold conditions expression of the four genes located in chromosome 3 was on average significantly repressed (Os03g07530, Os03g10750, and Os03g11420) or remained unchanged (Os03g12530) in lines with Guadiamar alleles and low V15 rates, while were on average upregulated (Os03g11420 and Os03g12530) or unchanged (Os03g07530 and Os03g10750) in cold tolerant RILs with the Italica Livorno haplotype ( Figure 6). It is interesting to note that in this second experiment gene expression under cold conditions was referred to that at 30 °C, the optimal for rice germination, instead of to gene expression at 25 °C. In this second experiment, differential expression of the Os11g05540 gene varied also between the two groups of lines, but patterns of expression were not consistent with those of the parent lines, particularly in lines with the Italica Livorno haplotype, in which expression under cold conditions was significantly repressed (mean Log2-fold change = −1.88 ± 0.79) while the gene was less repressed in Italica Livorno seedlings than in Guadiamar (Log2-fold change = −0.57 ± 0.03 front to −1.18 ± 0.07). This result was also inconsistent with that from the previous experiment (Figures 3 and 6), which can be related to the different temperatures used as control conditions. Therefore, our results indicate that the LTG phenotype in these japonica varieties would be explained by increased levels of expression under cold conditions of four genes (Os03g07530, Os03g10750, Os03g11420, and Os03g12530) from the rice chromosome 3.
One of these genes that in our two experiments showed contrasting differential expression patterns between cold-sensitive and tolerant rice lines has been previously associated with rice germination and cold stress. The locus Os03g07530 (gene Os03g0171600) encodes an F-box Kelch repeat domain-containing protein (OsFBK12) that was reported to interact with the S-adenosyl-L-methionine synthase 1 (SAMS1) and regulate ethylene biosynthesis [38]. In contrast with our results, overexpression of this gene in transgenic rice plants of the cultivar Zhonghua 10 delayed germination, while this process was accelerated in transgenic plants in which the OsFBK12 was repressed. F-box domain containing proteins are critical for the controlled degradation of cellular proteins by the ubiquitin/proteasome degradation pathway, where these proteins perform the crucial role of conferring specificity to the complex for appropriate targets [29]. Ubiquitination is a posttranslational modification affecting protein activity, localization, assembly, and interaction ability, that has been observed in many plant processes, including stress responses. The ubiquitin/proteasome system regulates phytohormone signaling, as occurs with ABA [39], since it has been referred to that of E3-ubiquitin ligases mediate ubiquitination of ABA receptors in several cell compartments. Recent works have described rice E3-ubiquitin ligases involved in plant response to low temperatures (see references in [40]), therefore this probably accounts for why several rice F-box proteins have been found to be differentially expressed under cold stress conditions [36,41]. Moreover, the candidate gene Os4g0619300 that encodes an F-box protein was mapped within a major QTL related to rice cold tolerance at booting [42]. From the genome sequence of Italica Livorno previously reported [43], we could determine that in this tagged region the Italica Livorno DNA sequence contained only one SNP, the selected marker in position 3,828,648, that is located 5291 bp upstream the gene, being the coding sequence of OsFBK12 that of the reference genome. A similar regulatory role is performed by the gene Os03g0205000, from the locus Os03g10750, that codifies a CUE (coupling of ubiquitin conjugation to ER degradation) domain-containing protein, Spotted Leaf 35 (SPL35), which is constitutively expressed in all organs and has a role on regulation of cell death and defense response via H2O2 accumulation [44]. The sequence of the Italica Livorno genome contains a base deletion 1400 bp upstream of this gene.
Significantly higher differential expression in rice cultivars and RILs showing the LTG phenotype was also observed for the gene Os3Bglu6 (Os03g0212800) that in locus Os03g11420 codifies for a chloroplastic β-glucosidase that regulates cellular ABA pools, affecting drought tolerance and photosynthetic activity [45]. Therefore, our results suggest that fine tuning of ABA levels would be the basis of sustained seedling growth during japonica rice germination, as was proposed in the previous studies [35]. The SNP marker from our panel (position 5,874,917) is 1327 bp upstream from the sequence of the Os03g11420 gene, and the Italica Livorno allele also differed from that of Nipponbare at the position 5,884,737 (1433 bp downstream the gene location), where an insertion was detected. We also found contrasting patterns of differential expression under cold conditions for the locus Os03g12530 (gene Os03g0226400) that encodes a Mn-transporter protein, OsMTP8.1, which is localized at the tonoplast [46]. Mn is essential for most photosynthetic organisms but its toxicity can also be detrimental to plants, therefore several mechanisms have arisen to deal with eventual excess of this micronutrient, such as Mn sequestering in the vacuole. The Italica Livorno sequence for Os03g12530 contains two SNPs in bases 2976 (within an intron) and 3318 (within the 7th exon but synonymous).
Germination under cold conditions increased the expression of genes Os03g12820 and Os06g06400 in the five cultivars studied (Figure 3), but Log2-fold change rates were higher in seedlings of LTG cultivars than in Guadiamar. However, expression of genes Os03g12820 and Os06g06400 was upregulated in RILs carrying the Italica Livorno allele for the tagged SNPs, while in RILs with Guadiamar alleles cold increased the expression of Os03g12820 at lower extent, and repressed the expression of Os06g06400 (Figure 7). Therefore, these two genes may have been also involved in rice LTG regulation. In fact, gene Os03g0230300 in locus Os03g12820 corresponds to the abiotic stress response OsSRO1c, an SRO homologue that regulates stomatal closure and abiotic stress tolerance by modulating oxidative stress [47]. This gene was induced by multiple stresses, including cold. When compared to Nipponbare, there are two mutations in the Italica Livorno allele of the gene 03g12820: an SNP 2893 bp upstream, and an insertion 705 bp downstream. The locus Os06g06400 (gene Os06g0158500) codifies an NBS-LRR type disease resistance protein that was reported to be overexpressed under abiotic stress [29], and was included in our study although the SNP in this location was not in the panel with stronger association to LTG.
The locus Os07g47470, that was also induced by cold in the five analyzed rice cultivars, encodes a pentatricopeptide repeat (PPR)-containing protein, while the PPR protein encoded by the gene Os02g0824000, from the locus Os02g57800, was upregulated at 15 °C only in seedlings from the LTG cultivars. In the experiment performed with RILs, Os02g57800 was repressed under cold conditions, particularly in Guadiamar and in RILs carrying the allele of this parent, while Os07g47470 was induced by cold in these lines and in Italica Livorno, and it was repressed in RILs carrying the Italica Livorno allele for this SNP and in Guadiamar ( Figure 6). Again, differences observed between experiments for cold-induced Log2-fold change values in parental lines could be due to the different temperatures used as reference (25 or 30 °C), but a consistent association between LTG phenotypes and gene expression under cold conditions was found only for the PPR-protein encoded by the Os02g0824000 gene. PPR proteins are involved in RNA processing, such as editing, splicing, stabilization, and translation, generally in chloroplast or mitochondria. During rice anaerobic germination, a number of genes encoding PPR domain containing proteins were found to be specifically induced [48], and also a PPR protein regulates growth and chloroplast development under cold stress [49]. The Italica Livorno alleles for Os02g57800 and Os07g47470 contain several mutations (SNPs and INDELs) in both upstream and downstream the coding sequences, while the allele for Os07g47470 contains also an SNP within the coding sequence (position 28,373,683) that causes the substitution of an asparagine residue by a threonine residue in the encoded protein. Similar results with inconsistent expression patterns in the two experiments performed were obtained when we analyzed those of a putative gene expression regulator, a GRAS family transcription factor (OsGRAS40) encoded by the gene Os10g0551200 from locus Os10g40390 (Figures 3 and 6), therefore the LTG phenotype could not be clearly associated with the cold induced differential expression of this gene. As mentioned earlier, this was the case also for gene Os11g0153400 in locus Os11g05540, that encodes a rho-GAP domain containing protein, the Rac GTPase activating protein 1. Rho-like GTPases from plants are involved in controlling cellular processes such as regulation of cytoskeleton and vesicular trafficking, and are positively regulated by auxins and negatively regulated by ABA [50].
The co-expression analyses (ATTEND-II) of the studied genes showed relationships with the metabolism of tryptophan, pyruvate and carbon, as well as with signaling pathways against biotic stresses (probable calcium-binding protein 11 and 48-CML11, CML48and MAPKKK), and mRNA regulation (LOC4334045 and LOC4343740). In this sense, a recent study identified several genes involved in the metabolism of tryptophan, a precursor of melatonin, in the association between the circadian rhythm and the response to cold stress in rice [51]. Taking these results together, we can conclude that LTG in these japonica rice cultivars was associated with transcriptional changes in five genes from chromosome 3 (Os03g07530, Os03g10750, Os03g11420, Os03g12530, and Os03g12820), as well as to cold-induced differential expression of genes Os02g57800 and Os06g06400. Therefore, the LTG phenotype was explained by the activation of genes that modulate ethylene and ABA levels (OsFBK12, Os3Bglu6), as well as regulatory genes (SPL35, Os02g0824000, and Os06g0158500) and genes involved in oxidative stress tolerance (OsSRO1c) and in Mn homeostasis maintenance (OsMTP8.1), rather than being explained by differences in the coding sequence of these genes. Enhanced gene expression rather than natural polymorphism in the coding sequence of OsbZIP23 explained drought tolerance in rice genotypes [52], while polymorphism within the coding sequence that implied changes at the protein level was detected in the candidate gene Os01g0620100, which showed different gene expression levels between cold tolerant and sensitive rice landraces germinating under cold stress [53].
Most of the LTG-associated SNPs are located in intergenic regions, as have been described for GWAS performed for other agronomic traits [54], and therefore might be part of the regulatory information contained within the non-protein coding parts of DNA, such as the long non-coding RNAs (lncRNAs) that are especially numerous in plants and a major regulator of the transcriptional process, including abiotic stress responses. For example, some rice lncRNAs have been found to be differentially regulated under drought and cadmium stress [55,56]. In the rice genome, intergenic lncRNAs are located close to flanking genes, and its expression correlated with that of the neighbor protein-coding gene in a genome-wide scale [57]. Two of the genes analyzed in this study (Os03g07530 and Os10g40390) are the closest ones to two lncRNAs detected in that work (LNC_Os03g07370 and LNC_Os10g02850, respectively). Alternative polyadenylation, which generates different transcripts with altered coding capacity, could also account for post-transcriptional regulation of LTG, since seven (Os02g57800, Os03g07530, Os03g11420, Os03g12530, Os07g47470, Os10g40390 and Os11g05540) out of the ten analyzed genes present more than a poly(A) site [58]. Finally, LTG might be regulated also by transposable elements, since these loci are transcriptionally activated in response to stress along with their neighboring genes [59]. Our study initially included other genes located in the regions from chromosomes 6 and 7 that were also associated with LTG by the previous GWAS [16], which encode transposon proteins, such as Os07g23290, which contains the SNP at position 13,128,278. Unfortunately, we could not obtain specific amplifications for these genes.
This complex regulation of LTG probably accounts for the low genetic improvement observed in V15 rates of individuals from the two progenies studied when inherited the haplotype of the donor parents, and particularly in RILs of the family Guadiamar × Italica Livorno (Figure 4b). In our experiments, Guadiamar showed low seedling growth rates at both suboptimal (V15 = 0.31 ± 0.07 mm/d) and optimal temperatures (V30 = 1.85 ± 0.28 mm/d), suggesting that improving LTG of this cultivar is limited by further intrinsic factors. However, substantial improvements in V30 rates could be observed in RILs ( Figure  4a) derived from the cross with Italica Livorno, which showed a higher growth rate under optimal temperature conditions (V30 = 3.36 ± 0.21 mm/d). Furthermore, in the Guadiamar × Italica Livorno F6 family, and in F3 lines derived from the cross Copsemar 7 × Nuovo Maratelli, significant increases in seedling growth rate under cold conditions were induced in plants carrying the donors' alleles for LTG-associated SNPs from chromosome 3 ( Figures 5 and 7), which tagged genes with differential expression patterns consistently associated with phenotypic variation for this trait. Therefore, the observed LTG genetic gain in these biparental families was mainly explained by introgression of favorable alleles in genes regulating ethylene and ABA metabolism and oxidative stress. High heritability estimates for LTG were reported [60], with both additive and non-additive gene actions involved in the regulation of this trait, and with the non-additive component relatively more important for percentage of reduction in coleoptile length and coleoptile growth. Epistatic interaction (a non-additive effect) was also found to be important for rice germination capacity under low temperature [61]. The importance of relatively high non-additive effects suggests that selection should be applied in advanced generations of breeding programs (F4 or F5), when dominance effects have been decreased due to the increase in homozygosity of the individuals. Despite this limitation, phenotypic and marker-assisted selection in the families Guadiamar × Italica Livorno and Copsemar 7 × Nuovo Maratelli, allowed us to obtain advanced lines with more rapid germination at both standard and adverse temperature conditions. Promising plant materials were further subjected to field evaluation for grain yield, other agronomic traits and grain quality for possible release, while several advanced lines will be registered.

Conclusions
Gene expression and genotypic analyses performed along our breeding scheme indicate that low temperature germinability in japonica rice can be partially explained by the activation of the expression of genes modulating ethylene and ABA levels, as well as genes involved in the response to oxidative stress and in Mn homeostasis maintenance. Despite the complex genetic architecture of the trait, it was possible to obtain improved japonica rice lines that incorporated favorable alleles for previously described QTLs located in chromosome 3. These results contribute to enlighten rice LTG regulation and describe genetic resources and QTLs that can be useful for rice breeders in temperate areas.
Supplementary Materials: The following are available online at www.mdpi.com/article/10.3390/agronomy11112125/s1, Figure S1: Seedling growth rates (mm/d) of five japonica rice cultivars. Figure S2: Rice seedlings of six F6 lines Guadiamar × Italica Livorno germinated at 15 °C. Table S1: Primers used for real time amplification of ten candidate genes.

Data Availability Statement:
The data presented in this study are available within the article.