QTL mapping and candidate gene analysis of low temperature germination in rice (Oryza sativa L.) using a genome wide association study

Plant materials

In this study, a collection of 276 rice accessions were selected from the 3K-RG project. All rice accessions were cultivated in the same geographical location in Huai’an (119°0′14″E, 33°38′43″N), Jiangsu province in 2019 and 2020. Each accession was subject to the same field management in 2019 and 2020. To eliminate error results caused by marginal effects, every rice accession was planted in a 5 × 5 block within the 3 m × 3 m square, and five plants of each accession were randomly chosen from the middle of each square as the experimental subjects.

LTG measurement

The seeds of each rice accession were collected independently in a nylon bag with dense nets to air dry seeds for 2 weeks. After that, air-dried seeds were placed in the oven at 50 °C for 7 days to break primary dormancy. A total of 100 plumped seeds of each rice accession were extracted and spread on a round wet filter paper and kept at 15 °C and in darkness for germination. The number of germinated seeds was recorded daily from D6 to D10 with a seed shoot or root exceeding 0.1 cm considered a germinated seed (Wang et al., 2018c; Akhtamov et al., 2020; Najeeb et al., 2020). Seed germination rate = germinated seeds/100. LTG was assessed according to the germination rate of each recorded day.

SNP filter analysis

The genetic variations of 276 rice accessions are available publicly in the 3K-RG database and the information for all SNPs can be downloaded from the website for free (https://snp-seek.irri.org/_download.zul). In this study, the set criteria for selecting high-quality SNPs were based on (1) minor allele frequency (MAF) ≥0.05 and (2) number of accessions with minor alleles ≥6 (Yang et al., 2014). After filtering, only high-quality SNPs were retained. A slide window of 1 Mb was adopted to demonstrate the distribution of variants in all 12 chromosomes to determine the density of the SNPs. The detected SNPs were annotated and the possible effects were predicted through ANNOVAR (Wang, Li & Hakonarson, 2010).

Population structure analysis

To analyze the population structure, a principal component analysis (PCA), a neighbor-joining (NJ) tree and a K value analysis were applied. The phylogenetic tree was constructed using MEGA7 (version 7.0) (Kumar, Stecher & Tamura, 2016) and the results were visualized using ggtree (version 1.7.10) (Yu et al., 2017). The PCA was conducted by PLINK (version v1.90) (Purcell et al., 2007). According to the Bayesian Markov Chain Monte Carlo (MCMC) Program, the K value, ranging from 2 to 7 in the full population was inferred using STRUCTURE (version 2.3.4) (Pritchard, Stephens & Donnelly, 2000). The optimal K value was determined by ΔK (Evanno, Regnaut & Goudet, 2005). The result was visualized and the relevant Q matrix was generated for further analysis.

Programs for GWAS analysis

Based on the factored spectrally transformed linear mixed model, two programs, FaST-LMM (version 0.5.1) and GEMMA (version 0.98.1), adding different genetic similarities to analyze random effects, were applied to perform the GWAS. The validated number of SNP markers (N) was calculated using the Genetic type I Error Calculator (GEC) software (Turner, 2014) and suggestive (1/N) P value threshold was adopted as the standard to control type I error.

Quantitative real-time PCR assay

Ten seeds of each accession were sampled at 15 °C and in darkness. Total RNA was extracted using the TIANGEN RNAprep Pure kit (#DP441; TIANGEN, Beijing, China) according to the manufacturer’s protocol. Complementary DNA (cDNA) was synthesized using a cDNA synthesis kit (#RR047A; Takara, Tokyo, Japan). Quantitative real-time PCR (qRT-PCR) reaction was conducted using TB Green Premix Ex Taq (#RR820A; Takara, Tokyo, Japan). The reaction was performed on the CFX Connected Real Time System (Bio-Rad, Hercules, CA, USA). The expression level was calculated by 2^−ΔCt using the expression level of Ubiquitin as reference. Each sample was tested three times to fulfill technical replications. Relevant primer sequences are provided in the supplemental data (Table S1).

Phenotypic variation for LTG in natural rice accessions

A collection of 276 rice germplasms was selected from the 3K-RG project for the LTG test. Rice accessions in this study were from 17 different regions worldwide (Fig. S1, Table S2). Previous studies have applied different temperatures ranging from 12 °C to 15 °C to estimate LTG (Borjas, De Leon & Subudhi, 2015; Fujino et al., 2008; Li et al., 2013; Wang et al., 2011). Given the effect of secondary dormancy induced in 12 °C (Miura & Araki, 1996), 15 °C was applied to evaluate LTG in this study. Germination was defined as the seed shoot or root exceeding 0.1 cm from the seed coat (Fig. 1A), and the evaluation of LTG was based on the germination rate from D6 to D10. The average germination rate on 2019D6 was lower (3.5%) than that (13.6%) on 2020D6 (Fig. 1B, Table 1), suggesting that environmental factors have an impact on the phenotype of LTG. Furthermore, the total number of germinated accessions was too small to draw any general conclusions (Fig. S2). The average germination rates were 22% and 31% on 2019D7 and 2020D7, respectively (Fig. 1B, Table S3). On 2019D7, a total of 129 rice accessions had a germination rate that exceeded 10% while on 2020D7 this number rose to 161 (Fig. 1C). From D8 to D10 in both years, the average germination rate almost reached or exceeded 40% (Fig. 1B, Table 1). The investigation of germination rate gaps between two adjacent days indicated that the gap between D6 and D7 in both years was the largest at 18.5% in 2019 and 17.4% in 2020, respectively (Table 1). In contrast, the smallest germination gap in 2019 was 13.8% between D8 and D9, while in 2020, the gap between D9 and D10 was the smallest (10.3%) (Table 1). On D10 of both years, all rice accessions had a relatively high germination rate ranging from 64.7% to 69.4%. Overall, the distribution of germination rate followed similar trends in both years although the germination rate of D6 and D7 in 2020 was higher than that in 2019 (Fig. 1B, Table 1).

SNP density analysis

The original version of the 3K-RG database contained 32 million SNPs in total. Through filtering, a total of 497,231 SNPs were detected. After the classification of SNPs, the density of SNPs in all 12 chromosomes were between 1,033.3/1 Mb and 1,648.94/1 Mb (Fig. S3, Table S4). This indicated that the filtered SNPs in this study were sufficient and distributed evenly in 12 chromosomes.

Population structure and kinship

Using the SNPs, we performed a PCA to quantify the population structure of these 276 accessions. The total variance explained by PC1 and PC2 was 35.60% and 16.79%, respectively (Fig. 2A). Based on the Nei’s genetic distance (Nei, 1972), the NJ tree was plotted separating the full group into two groups (Fig. 2B). Meanwhile, using STRUCTURE, the peak of ΔK appeared when K = 2, suggesting that the full population could be divided into two subgroups (Figs. 2C and 2D). These two subgroups corresponded to Japonica and Indica (Table S2), which is consistent with the findings of Wang et al. (2018a).

GWAS for LTG in rice

A total of 136,276 validated SNPs (MAF ≥ 0.05) were used for the GWAS through the FaST-LMM and GEMMA models. The GEC was used to calculate the indicator P value, which gave 7.41E−6 as the suggestive P value (−log (p value) = 5.13). According to a previous study, the distance of two adjacent lead SNPs within 200 kb was considered one QTL (Lv et al., 2016). A total of 37 QTLs with 54 SNPs were found using FaST-LMM for LTG from D6 to D10 in both years whereas 107 QTLs with 159 SNPs were detected using GEMMA (Table 2, Tables S5–S7). Nearly half of the QTLs identified in FaST-LMM (15/37) were also identified in GEMMA, suggesting FaST-LMM had stricter criteria in controlling false positive association (Table S6).

Using FaST-LMM, there were 26, 8, 5, 4 and 4 QTLs detected from D6 to D10, respectively, in 2019 and 2020 (Table 2). There were seven QTLs characterized repeatedly in the total (Table 2). Compared with QTLs reported before, 13 QTLs were co-localized within the interval of 1 Mb in this study (Table 2). Among these QTLs, 12 QTLs were associated with LTG, five QTLs were associated with cold tolerance at the seedling stage, and four QTLs were associated with both LTG and cold tolerance at the seedling stage (Table 2). These results confirmed that the GWAS results in this study were reliable for further candidate analysis. The remaining 24 QTLs that had been uncharacterized before were considered novel QTLs for LTG. Among the novel QTLs, it was notable that qLTG1-2 was repeatedly detected on D7 in both years using FaST-LMM and GEMMA (Figs. 3A–3C). Moreover, this QTL was detected on 2020D6 as well using FaST-LMM and GEMMA (Fig. S4). Therefore, further analysis was focused on qLTG1-2 with this repeated lead SNP (Chr.1_13340259).

Haplotype and expression analysis of the candidate genes

To further locate the candidate genes of qLTG1-2, two public germination-related transcriptome datasets (SRP277875, GSE115371) were adopted. Dataset SRP277875 contained the expression data at different time points of two rice accessions for germination (Yang et al., 2020a) and Dataset GSE115371 provided the expression data of one rice accession under aerobic conditions for germination at various time points (Narsai et al., 2017). According to the published transcriptome data, 13 expressed genes located in the qLTG1-2 interval were found, which were then used for further comparison analysis of the expression levels in low temperature between high and low germination accessions (Fig. S5). LOC_Os01g23600, LOC_Os01g23705 and LOC_Os01g23850 failed to be amplified, suggesting they exhibit very low expression levels in the seeds. Eight genes, including LOC_Os01g23590, LOC_Os01g23610, LOC_Os01g23630, LOC_Os01g23640, LOC_Os01g23680, LOC_Os01g23710, LOC_Os01g23740 and LOC_Os01g23870, did not show obvious differences between high and low germination accessions (Fig. S6).

LOC_Os01g23580 was located 90 kb from the lead SNP and was associated with abiotic stress in a GO analysis. Furthermore, the homolog of LOC_Os01g23580 in Arabidopsis has been shown to be involved in the regulation of auxin transport and response to several abiotic stresses (Li et al., 2005; Wijewardene et al., 2020). One non-synonymous SNP (Chr.1_13243045, base G-C, amino acid Ser-Thr) and one upstream SNP (Chr.1_13236390, base A-G) were identified within the sequence of LOC_Os01g23580, which generated three haplotypes in the full population (Fig. 4A). Haplotype I of LOC_Os01g23580 displayed better performance for LTG than the left two haplotypes (Figs. 4A and 4B). Moreover, accessions of high germination rates were usually ones with G allele whose transcriptional levels were much lower than accessions of low germination rates with an A allele in the upstream region (Fig. 4C). LOC_Os01g23620, namely OsSar1a, was located 50 kb from the lead SNP and OsSar1abc RNAi mutants led to pre-harvest sprouting (Tian et al., 2013). Based on one upstream SNP (Chr.1_13285882 base A-G), the full population was divided into two haplotypes (Fig. 4D). Haplotype I of OsSar1a showed a higher germination rate than haplotype II which was negatively associated with transcriptional level (Figs. 4E and 4F).