Chilling tolerance trait phenotyping of RDP1 accessions
To better understand the natural variation that allows the JAPONICA subspecies of rice to be more chilling tolerant than the INDICA subspecies, multiple phenotypic analyses on a subset of Rice Diversity Panel 1 (RDP1; Eizenga et al. 2014) accessions were done at five different temperatures (4°C, 8°C, 10°C, 12°C, and 16°C). The core subset of 354 out of 424 RDP1 accessions was chosen due to their consistently robust germination rates and abundant seed availability, and because their genotypes were included in the High-Density Rice Array (HDRA) containing 700,000 single nucleotide polymorphisms (SNPs; Chen et al. 2014). A principal component analysis using the SNP information showed that there was indeed considerable genetic variation that led to a clustering of the different subpopulations, which is a critical requirement for genome wide association studies (GWAS) and quantitative trait loci (QTL) mapping (Fig. 1). The rationale for choosing different chilling temperatures was to determine whether different stress tolerance pathways are active at different temperatures, or whether common pathways can be identified that are active at overlapping temperatures. The low-temperature seedling survivability (LTSS) phenotypic assay was chosen, because it was previously shown that at 10°C, it significantly discriminated between accessions with different cold stress tolerance potentials (Schläppi et al. 2017). The electrolyte leakage (EL) phenotypic assay was chosen as a measure of plasma membrane integrity after cold temperature stress exposure, and it was previously shown that it can identify accessions with different cold stress tolerance potentials even when no differences in LTSS scores are observable (Shi et al. 2020). Phenotypic assays for EL and LTSS are shown in Fig. 2. Lastly, the median lethal chilling temperature (LT50), that is, the lethal temperature at which 50% of two-week-old seedlings die and 50% live after one week of cold exposure was calculated from LTSS scores to potentially identify additional chilling stress tolerance pathways at the threshold of survivability.
Collectively, as shown in Fig. 3 for three temperatures (8°C, 10°C, and 12°C), the JAPONICA and INDICA subspecies had significantly different LTSS, EL, and LT50 scores. Accessions within a particular subpopulation (aromatic, aus, indica, temperate japonica, tropical japonica, admixed JAPONICA, admixed INDICA, and admixed INDICA/JAPONICA) had similar chilling tolerance phenotypes, and they followed similar trajectories at every temperature tested. Interestingly, while most members of the JAPONICA subspecies were more chilling tolerant than members of the INDICA subspecies, aromatic accessions which are typically classified as part of the JAPONICA subspecies had chilling tolerance phenotypes that were more similar to admixed INDICA/JAPONICA than to other JAPONICA relatives, most likely because genotypically they also clustered with admixed varieties (Fig. 1). In agreement with our results, it was previously observed that aromatic accessions clustered with admixed instead of temperate or tropical japonica accessions (Wang et al., 2018). For this reason, aromatic and admixed INDICA/JAPONICA accessions were grouped into a third “Intermediate Cold Tolerant” (Intermediate) cluster, distinct from the “Cold Tolerant” (Tolerant; temperate japonica, tropical japonica, and admixed JAPONICA) and “Cold Sensitive” (Sensitive; aus, indica, and admixed INDICA) clusters. Admixed JAPONICA and admixed INDICA accessions predictably had phenotypes in accordance with their JAPONICA and INDICA mixed genotypes, respectively. However, it is interesting to note that admixed JAPONICA accessions were generally a bit less chilling tolerant than either temperate or tropical japonica accession while admixed INDICA accessions were generally a bit more chilling sensitive than either aus or indica accessions (Fig. 3). The three “chilling tolerance clusters” were used for multiple chilling tolerance trait analyses, which was the basis for subsequent GWAS mapping.
Summary of percent electrolyte leakage (EL) results
The percent EL was determined for each cultivar as a measure of plasma membrane integrity after cold temperature stress exposure, with the general expectation that cold sensitive accessions will release more electrolytes (high EL scores) from leaf tissues than cold tolerant accession (low EL scores), due to a higher incidence of membrane lesions after cold exposure. The mean percent EL for Tolerant, Intermediate, and Sensitive clusters generally followed this expectation (Fig. 4). For the Tolerant cluster, the mean percent EL at the five temperatures ranged from 8.7 to 20.9%; for the Intermediate cluster, 9.4 to 34.3%; and for the Sensitive cluster, 12.8 to 35.4%. At each temperature, mean EL values were significantly lower for the Tolerant cluster than for the Sensitive clusters. Mean EL values were also significantly lower for the Tolerant cluster than the Intermediate cluster at 4°C and 8°C, however, there were no significant differences between 10-16°C. This indicates that below 10°C, the Intermediate cluster shifts from a more tolerant to a more sensitive EL phenotype. For the Sensitive cluster, similar EL values at 4°C and 8°C suggests that 8°C is a threshold temperature for Sensitive accessions below which lower temperatures do not significantly increase membrane lesions any more. Taken together, these EL results suggest that Tolerant, Intermediate, and Sensitive clusters have different low-temperature thresholds at which maximum membrane lesions occur, and that each cluster has a different low-temperature buffering range for maintenance of membrane integrity.
Summary of low-temperature seedling survivability (LTSS) results
The percent LTSS was determined as a measure of survival of two-week-old seedlings after one week of continuous low-temperature stress exposure. The phenotype of plants after a 7-day recovery period from cold stress is show in Fig. 2: recovered seedlings were alive, green, and healthy looking while dead seedlings were bleached and wilted. Our expectation was that the Tolerant cluster will have significantly more accessions with high LTSS scores than the Intermediate and Sensitive clusters. The mean percent LTSS for Tolerant, Intermediate, and Sensitive clusters generally followed this expectation (Fig. 5). For the Tolerant cluster, the mean percent LTSS at the five temperatures ranged from 20.1 to 98.2%; for the Intermediate cluster, 3.0 to 96.2%; and for the Sensitive cluster, 0.3 to 94.9%. At each temperature, mean LTSS scores were significantly lower for the Sensitive cluster than for the other two clusters, and except for 16°C, mean LTSS values were significantly higher for the Tolerant cluster than the Intermediate cluster. Between 4-10°C, mean LTSS values for the Sensitive cluster were very similar, indicating that 10°C might be a threshold temperature for Sensitive accessions below which lower temperatures do not significantly decrease seedling survivability, which is similar to what was observed for the EL phenotype. Between 8-10°C, LTSS values for the Intermediate cluster shifted from a more tolerant to a more sensitive phenotype, which is also similar to what was observed for the EL phenotype, however, an apparent threshold for survivability was reached between 4-8°C. The somewhat different temperature threshold buffering ranges for mean EL and LTSS scores suggests that there may not be a linear correlation between membrane integrity and seedling survivability.
Correlation between LTSS and EL at different low temperatures
Percent EL and LTSS phenotypes of the Intermediate cluster were similar to those of the Tolerant and Sensitive clusters at high and low temperatures, respectively, but at different threshold temperatures (Fig. 4; Fig. 5). This suggests that the EL phenotype, measured immediately after one week of low-temperature stress exposure, might not accurately predict survival rates one week after recovery at growth promoting temperatures. Therefore, correlations between EL and LTSS were done to determine how effectively EL scores can be used to predict LTSS. Interestingly, linear correlation coefficients varied between the five temperatures (Fig. 6). The best correlation was observed at 8°C (R2 = 0.4620), followed by 10°C (R2= 0.3696) and 12°C (R2= 0.2869), while there was no correlation at the extreme ends of 4°C (R2= 0.0746) and 16°C (R2= 0.0098). This indicates that at 8°C ± 1°C, the degree of chilling temperature induced membrane damage correlates well with subsequent probability for survival. That is, less membrane damage predicts higher survival rates, however, even at 8°C, there were several outliers with little membrane damage and low survival rates, and several outliers with a high degree of membrane damage and high survival rates (Fig. 6). This suggests that for some accessions, membrane damage sustained during cold exposure is not the main cause for low survival rates, indicating that physiological or metabolic events during the recovery phase affect survival rates. In contrast, even though significant membrane damage occurs during cold exposure, some accessions can deal with this damage during the recovery phase, leading to increased survival rates.
Summary of median lethal chilling temperature (LT50) results
Our data show that different “cold tolerance clusters” (Fig. 4; Fig. 5) as well as different subgroups (Fig. 3) had different “tipping points” where EL or LTSS scores changed dramatically. To include these tipping points as an additional parameter for GWAS analyses, we calculated LT50 values from LTSS data for each accession at which 50% of seedlings died and 50% survived. Using R, LTSS curves were generated by plotting median LTSS rates at 4°C, 8°C, 10°C, 12°C, and 16°C to identify LT50 for each accession (Fig. 7a). The curves were clustered for Tolerant, Intermediate, and Sensitive accessions, and median LT50 were calculated for each cluster (Fig. 7b). The survival curves and median LT50 values for each tolerance cluster generally followed the patterns observed for LTSS, further supporting our classification of aromatic accessions as an Intermediate cluster.
Genome-wide association study (GWAS)-based mapping of EL, LTSS, and LT50 QTL
From the core subset of 354 RDP1 accessions, EL, LTSS, and LT50 data were used to perform GWAS to identify significantly trait associated SNPs. We applied a publicly available GWAS pipeline (McCouch et al. 2016) that accounted for population structure of the rice accessions using a kinship matrix created by the efficient mixed-model association eXpedited (EMMAX) approach (Kang et al. 2010), which determined ancestry or cryptic relatedness between individual genotyping markers. GWAS-mapping pipeline produced quantile-quantile (Q-Q) plots and Manhattan plots for each cold tolerance trait are shown in Fig. S1 and Fig. 8, respectively. Although several Q-Q plots appeared to indicate inflation in some of the phenotypes, our plots had similar trends as those of the original mapping results reported in the study that established the RDP1 mapping pipeline (McCouch et al. 2016). Therefore, based on Q-Q plots a p-value cutoff for significant SNPs was determined to be at -log10(p) > 4, and QTL were identified as chromosomal sites containing 3 significant SNPs within a one million base pair (Mbp) region.
Across the three traits and five temperatures, 245 total QTL were identified. The LTSS trait yielded 150 QTL, the EL trait 92 QTL, and the LT50 trait 3 QTL (Table S1-S7). Temperature also affected the number of QTL identified. The moderate chilling temperatures of 10°C and 12°C yielded 83 and 62 QTL covering 11.6% and 8.9% of the rice genome, respectively, while 4°C, 8°C, and 16°C yielded 38, 29, and 30 QTL respectively (Table S1-S5). Many of the 245 QTL were found in overlapping chromosomal regions, some of which were longer than 1 Mbp (Fig. 9; Fig. S2). Taken together, we identified 178 unique QTL, consisting of 40 multiple-trait QTL (qMT) and 138 single-trait QTL (qE for qEL, and qL for qLTSS) (Table S7), covering approximately 25% of the rice genome.
Using different cold tolerance clusters, additional GWAS-based QTL mappings were done. To identify QTL that control chilling tolerance in accessions of the Tolerant and Sensitive clusters, the Intermediate cluster was removed and GWAS mappings were done for the three high EL-LTSS-correlation temperatures of 8°C, 10°C, and 12°C (Table S8). GWAS mapping was also done within the Tolerant and Sensitive clusters to identify QTL correlating with phenotypic variability within each subspecies for the three high EL-LTSS-correlation temperatures of 8°C, 10°C, and 12°C (Table S9). The Tolerant versus Sensitive (TvS) GWAS mapping resulted in 219 individual QTL: 68 EL QTL, 150 LTSS QTL, and 1 LT50 QTL. The within-subspecies GWAS mapping resulted in 36 individual QTL: Mapping within the Sensitive cluster yielded 33 QTL consisting of 3 EL QTL, 29 LTSS QTL, and 1 LT50 QTL; and mapping within the Tolerant cluster yielded 3 QTL consisting of 2 EL QTL and 1 LTSS QTL. This result is consistent with the larger genetic variation within Sensitive (Fig. 1, red cluster) than Tolerant (Fig. 1, green cluster) accessions.
Within the 178 unique QTL regions obtained using all 354 accessions, 14,710 genes were identified, which is more than 25% of the total number of annotated genes reported by the Rice Genome Annotation Project (Kawahara et al. 2013). Since a major focus of our GWAS analysis was to identify general cold tolerance genes with natural variation between chilling tolerant and chilling sensitive rice accessions and to reduce the list of potential candidates, we interrogated the 40 qMT QTL regions for chilling tolerance candidate genes (Table 1). There were 30 QTL that overlapped with previously published cold tolerance QTL (Andaya and Mackill 2003; Liu et al. 2013; Ma et al. 2015; Wang et al. 2016; Lv et al. 2016; Schläppi et al. 2017). Genes found in the qMT regions are hypothesized to be involved in either multiple traits (EL, LTSS, LT50) or in helping with improving those traits at more than one temperature. This reduced the total number of genes to be analyzed from 14,710 to 6,065 genes associated with the 40 qMT regions, which were further analyzed using gene function and Gene Ontology (GO) term annotations.
To identify cold tolerance candidate genes within the qMT QTL, annotated genes that contained significant SNPs (p<10-4) within their genomic regions were selected, yielding a list of 220 genes. Significant SNPs with expected percentages of major alleles and direction of allele effect on chilling tolerance trait were used to further filter the gene list. That is, significant SNPs with a major allele frequency at least as much as the percentage of accessions of the Tolerant cluster (205 out of 354 accessions = 58%) and with an allele effect positively correlating with the measured chilling tolerance trait (to decrease EL and to increase LTSS) were selected, yielding a list of 81 candidate genes. Those genes were then used for GO term enrichment analysis and compared to the larger gene list to narrow down possible cold tolerance response pathways. Furthermore, since GWAS-based QTL mapping was also done using accessions from between and within the major cold Tolerant and Sensitive clusters, 35 qMT QTL regions with at least one single QTL using these clusters were identified and used for targeted GO enrichment analysis (Table 1). In addition to using a kinship matrix, this analysis was done to reduce the potential influence of population structure in defining significant regions. Thus, the rationale for using regions that were highlighted in mappings within clusters was to account for a potential population structure bias arising from only the differences in variability between the Sensitive and Tolerant clusters. Those 35 qMT QTL were further filtered for genes containing significant SNPs, which resulted in 207 genes that were separately used for GO term enrichment analysis.
GO term analyses of all genes within multi-trait (qMT) QTL
The qMT QTL gene list was purged of annotated transposable elements resulting in 4,256 genes which were analyzed for GO term enrichment. GO term annotations were obtained for 2,483 of the 4,256 genes, and 257 GO terms were found to be enriched (p<0.01), and those relevant to plant chilling tolerance were grouped into chilling stress relevant clusters. Within the GO term for Biological Process, three major clusters were identified for “carbohydrate derivative biosynthesis”, “protein modification by small protein conjugation or removal”, “carbohydrate transmembrane transport”, and “hydrogen peroxide mediated signaling pathway” (Fig. S3). Within the GO term for Cell Component, there were only two clusters: “transferase complex” and “intrinsic component of plasma membrane” (Fig. S4). Within the GO term for Molecular Function, numerous clusters were identified including: “ADP binding”, “carbohydrate transmembrane transporter activity”, “oxidoreductase activity”, “thiol-dependent ubiquitinyl hydrolase activity”, and “ubiquitin-like protein transferase activity” (Fig. S5). These terms suggest possible associations with chilling tolerance mechanisms and can be used to select pools of genes for further analysis. These genes are expected to be involved in metabolic homeostasis, membrane composition, ubiquitination, oxidative stress, and signaling pathways.
GO term analysis of filtered gene list from qMT QTL
The filtered list of 81 genes containing significant SNPs was also purged of annotated transposable elements, reducing it to 71 genes for GO term enrichment analysis. Due to the low gene number, significance of enrichment was set to p<0.05. GO term annotations were obtained for 45 genes and 10 terms were enriched. In Biological Process, “cellular response to stress” and “transmembrane transport” were enriched (Fig. S6). In Cell Component, “Golgi membrane” and “extracellular region” were enriched (Fig. S7). In Molecular Function, “molecular transducer activity” and “exonuclease activity” were enriched (Fig. S8). These enrichment clusters suggest likely interaction of genes involved in vesicle transport and secretion and nucleic acid processing to confer cold stress tolerance responses. Enrichment for cellular response to stress and Golgi membrane suggest that these filtered genes may be involved in modifications and transport of membranes and membrane-bound proteins.
GO term analyses of significant genes within multiple-trait and cold Tolerant & cold Sensitive cluster QTL
Genes within the 35 qMT QTL that overlapped with cold Tolerant and cold Sensitive cluster-specific QTL (Table 1) were separately investigated to identify enriched GO terms. Out of 67 input genes, the full filtering protocol resulted in only three enriched (p<0.05) GO terms: “exonuclease activity”, “cytoplasmic part”, and “extracellular region”. As term enrichment map visualizations were not possible with so few terms, an assessment was done on a larger list of genes from an intermediate filtering protocol. The intermediate filtering of only genes with significant SNPs resulted in 159 non-transposable element genes with significant SNPs from the 35 qMT and Tolerant & Sensitive Cluster QTL. The terms that were enriched at this stage were related to “ubiquitination”, “plant growth”, and “nucleotide binding” (Fig. S9-S11).