Whole-genome meta-analysis coupled with haplotype analysis reveal new genes and functional haplotypes conferring pre-harvest sprouting in rice

Abstract

Background

Pre-harvest sprouting (PHS), which adversely impacts grain yield and quality, is controlled by seed dormancy genes. However, only a few dormancy-related genes have been characterized, and the effects of allelic variation in genes and the genetic basis of seed dormancy in rice remain largely unknown. Here, we performed a whole-genome meta-quantitative trait loci study to elucidate the genetic basis of seed dormancy in rice.

Result

One hundred and sixty-seven QTL were identified for PHS from which 134 were successfully projected onto the reference map yielding 20 consensus regions, meta-QTL (mQTL). The mean confidence interval of the mQTL was narrower (9.56-fold reduction) than that of the initial QTL. Six of the 20 identified mQTL were designated as breeders’ mQTL based on their small confidence intervals, large phenotypic variance explained, and the involvement of high number of QTL. Further, we retrieved 559 high-confidence genes from breeders’ mQTL regions conferring resistance to PHS. Comparative analysis of genes found in breeders’ mQTL loci and an RNA-seq-based transcriptomic dataset discovered 34 common genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed a significant enrichment of the common genes in amino sugar and nucleotide sugar metabolism, carbon metabolism, and carbon fixation in photosynthetic organs. Combined in silico expression profiling and qRT-PCR validation showed that LOC_Os10g18364, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140 exhibited high fold-change expression in PHS resistant cultivar (23xS-261) than PHS susceptible cultivar (23xS-262). Association analysis of these genes with germination rate index demonstrated that LOC_Os10g18364^Hap1, LOC_Os10g21940^Hap1, LOC_Os10g22590^Hap1, and LOC_Os10g25140^Hap1/Hap3 exhibited low germination rate (GR) in cultivars carrying these haplotypes.

Conclusion

In summary, this study delineates the genetic basis of PHS and provides a new set of target genes for improving PHS resistance. The natural variants identified in these genes and markers associated with breeders’ mQTL serve as potential resources for incorporating PHS resistance in rice.

Background

Dormancy in cereal seeds is a survival strategy that regulates the timing of germination under intricate environmental conditions. However, inadequate seed dormancy can trigger a phenomenon referred to as pre-harvest sprouting (PHS). Crops that exhibit reduced dormancy have been selected by humans during domestication, which allows for improved crop management, but at the same time increases the risk of seed sprouting before harvest [1,2,3]. Under high temperature and humidity, PHS in rice reduces yield and seed quality posing financial setbacks for farmers and rice industries across the globe [4]. Also, PHS affects cooking, processing and eating qualities due to changes in grain traits, starch structure and physicochemical properties [5, 6]. A deteriorating climate further poses more challenges to rice cultivation, thus breeding PHS resistance in rice cultivars is imperative.

Seed dormancy is a complex quantitative trait governed by numerous quantitative trait loci (QTL)/genes and influenced by internal and external environmental cues [7]. Phytohormones-mediated signaling pathways are involved in PHS and seed dormancy, such as abscisic acid (ABA), gibberellic acid (GA) and ethylene [8, 9]. Key regulatory genes including Sdr4, OsDSG1, OsAB13, OsAB15, PHS8, PHS9, OsNCED3, OsVP1, OsPDS, β-OsLCY, OsFbx352, OsMFT2, OsZDS, and OsCRTISO are known to modulate seed dormancy and germination via ABA synthesis, perception, catabolism and signal transduction [10,11,12,13]. Recently, OsGA2ox9 a gene encoding a GA 2-oxidase was shown to promote seed dormancy by altering GA metabolism and reducing ABA signaling [14]. To date, several QTL governing rice seed dormancy have been reported [15,16,17], however, only a handful of these QTL have been position-cloned and our knowledge about the genetic and physiological factors underlying PHS resistance remains limited. This calls for an augmented effort toward dissecting genomic regions/QTL and identifying more regulatory genes responsible for PHS.

The mapping of QTL underlying seed dormancy in segregating populations from different genetic backgrounds has unraveled the genetic basis and detected many QTLs [18, 19]. However, variations in factors such as experimental design, experimental environment, genetic background of parental materials, population type and size, molecular marker density, and statistical methods employed generally yield varied QTL positions which are inconsistent and often display low phenotypic variance. Thus, further analysis is required to clarify the ambiguity associated with identified QTL and pinpoint their location precisely. The meta-analysis is a reliable approach that provides more robust QTL including identifying QTL hotspots in a large dataset [20]. The meta-QTL (mQTL) approach has been successfully employed to identify mQTL in a range of traits including grain yield and related traits [21,22,23,24], biotic and abiotic stress tolerance [25,26,27], nitrogen use efficiency and root system architecture development [28, 29], and quality traits [30, 31].

The pioneering work of Tyagi and Gupta [32], on mQTL for PHS in wheat compiled QTL results from 15 studies involving 15 diverse mapping populations and identified eight meta-QTL with reduced confidence intervals (CI). Surprisingly, the mQTL positions overlapped with wheat viviparous 1 (taVP1) and Giberrellin 20 oxidase 1-A (TaGA20-ox1) genes which play a vital role in regulating PHS resistance in wheat [32,33,34,35]. In another study, Tai et al. [36] integrated 188 original QTL from 40 different QTL studies and confined them into 66 mQTL distributed on all 21 wheat chromosomes which revealed signaling networks and environmental factors that underpin PHS/seed dormancy. However, no meta-QTL study has been conducted on PHS/seed dormancy in rice. To the best of our knowledge, this is the first-ever attempt to employ mQTL analysis approach to identify mQTL regions associated with traits contributing to PHS resistance in rice.

In this study, we exploited existing QTL data published in previous literature to unravel the genetic and physiological factors influencing PHS resistance in rice. In particular, the objectives were: (i) to conduct mQTL analysis of QTL underlying PHS resistance in rice, (ii) to mine candidate genes within mQTL loci and identify novel genes regulating PHS, (iii) to perform association analysis to identify functional markers associated with PHS resistance in rice. Our study is yet another effort that provides baseline information for understanding the genetic and physiological processes toward genetic improvement of pre-harvest sprouting resistance and or seed dormancy in rice.

Materials and methods

Literature search, QTL data collection, and input file Preparation

A thorough bibliographical search for PHS/seed dormancy-associated traits QTL published in a span of 23 years (1998 to 2023) was carried out on Google Scholar (https://scholar.google.com/), PubMed (http://www.ncbi.nlm.nih.gov/pubmed), and gramene database (https://archive.gramene.org/qtl/) using appropriate keywords. The salient features of these QTL are summarized in Table 1. For each QTL study, the collected QTL information included: (i) QTL names/ID; names of associated traits; (ii) logarithm of odd (LOD) values of each QTL; (iii) phenotypic variance explained (PVE) or R² values; (iv) flank markers positions; (v) type and size of mapping populations (F₂, Backcross population, double-haploid population, and Recombinant inbred lines). QTL with incomplete information were ruled out from further analysis to guarantee data integrity. In instances where PVE and LOD scores of QTL were missing, 3% and 10% were used based on the proponent of Venske et al. [37]. Also, the 95% confidence interval of each QTL was recalculated following the standardized formula: (1) F₂ and backcross population, CI = 530/(Mapping population × R² ); (2) RIL, CI = 163/(Mapping population x R² ); where R² denotes the phenotypic variance explained values of QTL and mapping population comprised the number of lines used in the initial QTL studies. The details of these QTL are presented in Additional file 1: Table S2.

Table 1 List of QTL mapping studies used in the meta-analysis

Construction of the high-resolution consensus map and QTL projection

For the construction of the integrated consensus linkage map, three reference maps from previous studies were used; (i) [38]; (ii) [25]; (iii) [39] (Fig. 1). Moreover, marker data from each study were also incorporated into the consensus map. The LPMerge package [40] in RStudio software was utilized to construct the reference map. A two-step approach was adopted as described [41]. The initial step is to determine the number of markers, the number of consensus bins, and the first count of ordinal conflicts. Next, the package provides statistical techniques for resolving conflicts when markers display inconsistent order. The second step involves generating a series of reference maps (K = 1 to 4), where K denotes the maximum interval size. In addition, the root-mean-square error (RMSE) between each map and the reference map was estimated. The reference map that exhibits the lowest RMSE was selected for subsequent analysis. For the QTL projection, two approaches were considered. First, the method proposed by Goffinet and Gerber [20] was used when the number of QTL on each chromosome is below 10. Following this, the best model with the lowest Akaike information criterion (AIC) score in BioMercator (https://versailles.inra.fr/Tools/BioMercator-v4) was selected. For the second approach, a two-step method proposed by Veyrieras et al. [42] was used when the number of QTL on individual chromosomes is above 10. The mQTL in each chromosome were then determined based on AIC model, corrected AIC (AICc and AIC3) model, Average Weight of Evidence (AWE) model, and Bayesian Information Criterion (BIC) model. The model with the lowest values was selected for subsequent analysis. The 95% CI and mQTL positions were determined following the best-selected model in the previous step. The QTL were then integrated to obtain the peak position of the initial QTL involved in mQTL CI. However, QTL that did not meet the minimum AIC threshold were discarded. mQTL were named based on their genetic positions on the rice chromosomes (i.e. mQTL1.1, mQTL2.1, mQTL3.1). The PVE of mQTL was calculated as the average PVE of the initial QTL involved in each mQTL.

Fig. 1

Venn diagram of unique and common markers among the three maps utilized for consensus map development

Validation of peak markers flanking breeder mQTL region

The peak Simple Sequence Repeat (SSR) marker RM496 and Restriction Fragment Length Polymorphisms (RFLP) marker RG207 associated with mQTL5.1 and mQTL10.1 were selected for validation on 16 rice genotypes. All genotypes were obtained from the seed bank of Yangtze University, College of Agriculture. Genomic DNA was isolated from fresh leaf samples of all genotypes using the standard Cetylmethyl ammonium bromide (CTAB) protocol with minor modification [43]. The following ingredients were utilized to produce CTAB: 2% CTAB, 1.5 M NaCl, 10 mM Na₃EDTA, 0.1 M HEPES-acid, 100% isopropanol, 70% ethanol, and 1xTE (10 mM of Tris-HCl, pH 8.0; 1 mM of EDTA). For both markers, PCR was performed in 20 µL using 10X PCR buffer, 25 mM MgCl₂, 10 mM dNTPs, 2 µM primers, 5 U Taq DNA polymerase, and 100 ng of template DNA. The following PCR cycling conditions were used: denaturation at 94 °C for 30 s, annealing at 60 °C for 30 s, extension at 72 °C for 10 s, and a final extension at 72 °C for 25 min. The amplified products were separated by electrophoresis in a 1.2% agarose gel and visualized using the Gel Documentation System (BIO-RAD 2000).

Survey of transcriptome to identify candidate genes within breeders’ mQTL intervals

The flanking markers of breeders’ mQTL were manually searched and primer sequences were retrieved from the Gramene Marker Database (http://archive.gramene.org/markers). The primer sequences of the markers were subjected to BLASTN against the rice reference genome in phytozome (https://phytozome-next.jgi.doe.gov/) to retrieve the physical positions. For markers without physical coordinates, the physical positions were manually anchored. Subsequently, genes in these mQTL regions were batch-downloaded from Rice Genome Annotation Project Database using the physical position of the markers (http://rice.uga.edu/cgi-bin/gbrowse/rice/). A thorough survey of differentially expressed genes (DEGs) from RNA-seq study (Accession ID: PRJNA1020757) for seed dormancy in rice [44] was used to obtain common genes between DEGs and candidate genes (CGs) in breeders’ mQTL regions. The criteria employed for DEGs between genotypes were a log2-fold greater than 1 and a p-value less than 0.05 (p < 0.05). The Venny online tool (https://bioinformatics.psb.ugent.be/webtools/Venn/) was used to compare the DEGs and the genes in breeders’ mQTL regions. To further comprehend the functional attributes of the common genes, Gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) were performed using Comprehensive Annotation of Rice Multi-omics tool (http://bioinfo.sibs.ac.cn/carmo/Gene_Annotation.php) and the KOBAS v3.0 software (http://kobas.cbi.pku.edu.cn/). Further, we visualized the results using TBtools software and bioinformatics online tool (http://www.bioinformatics.com.cn/), respectively.

In Silico expression profile of common candidate genes and protein-protein interaction network analysis of priori candidate genes

The In-silico expression of common genes was evaluated using a public transcriptome dataset available on the Rice Genome Annotation Project Database (http://rice.uga.edu/expression.shtml). The dataset consisted of differentially expressed genes in endosperm (25 days after pollination (DAP)), seed, and embryo (25 DAP). In this case, the expression data was presented as log2 transformed transcript per million (TPM) values. Only genes with substantial fold change (FC) ≥ 2 were considered DEGs. TBtools software [45] was used to visualize the gene expression patterns. Also, developmental atlas expression data from the Rice eFP browser (https://bar.utoronto.ca/eplant_rice/) was further used to explore the expression of candidate genes. The protein sequences of priori candidate genes (CGs) were downloaded from phytozome database and used to predict possible interacting proteins via STRING database (https://cn.string-db.org/) using the combined score > 0.9 and visualized using cytoscape software [46] (http://www.cytoscape.org/). The annotations of the priori CGs interacting proteins were obtained from gene ontology (GO) database (http://www.geneontology.org/).

Plant materials, growth conditions, and seed germination assay

Two rice cultivars that showed distinct polymorphism at the validated loci were selected for germination assay. The resistance and/or susceptibility of these two rice genotypes to PHS were re-validated. A total of 60 detached seeds harvested 45 days after heading (DAH) from each rice genotype were placed on moist paper in Petri dishes and kept at 100% relative humidity (RH), and 27 °C in a growth chamber regulated at 16 h light/8 h dark for 4 days (96 h in total). The germination rate (GR) of detached grains on the 4th day was calculated as follows: (number of seeds germinated/total filled seeds) × 100% [47]. In line with the observed polymorphism, the two cultivars showed significant differences in germination rate and therefore were used for gene expression analysis of priori CGs for PHS resistance.

RNA extraction, cDNA synthesis, and expression analysis

Freshly collected embryos subjected to 12 h imbibition were snap-frozen in liquid nitrogen and kept at -80 °C for subsequent analysis. RNA extraction was performed using RNA prep Pure Plant Kit (TIANGEN, China). Total RNA was quantified using Nanodrop ND-25 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA). The RNA quality was thoroughly assessed through electrophoresis on a 1.2% denaturing agarose gel. First-strand cDNA of 1 µg total RNA was reverse transcribed using Takara PrimeScript Kit (Takara Bio). An SYBR Premix Ex TaqTM Kit was used to perform a quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysis on a Bio-Rad CFX96™ real-time system. The Rice Actin1 gene was used as internal control. The cycle threshold (Ct) comparative method (2^−ΔΔCt) was employed to quantify the relative expression level of priori genes, which is expressed as n-fold changes in transcript levels that were normalized with regard to an endogenous Actin gene [48]. The list of primers used is presented in Additional file 1: Table S6.

Haplotype analysis of CGs involved in PHS resistance

The SNPs in intragenic regions of the priori candidate genes resequenced from 50 selected rice cultivars consisting of 42 indica and 8 japonica were obtained from the 3 K Rice Genome Project from the International Rice Genebank Collection Information System (http://www.snp-seek.irri.org) (Additional file 1: Table S7). To verify their association with PHS resistance, the detached seeds of each rice accession harvested at 45 DAH were subjected to 100% relative humidity and 27 °C in a growth chamber at 16 h light/8 h darkness (96 h in total). The GR of each accession was recorded and used for association analysis (Additional file 1: Table S7). For haplotype analysis, all single nucleotide polymorphism (SNP) markers from the priori genes were utilized except missing and heterozygote data, then the student’s t-test was performed to ascertain the significant variations among the haplotypes. Linkage disequilibrium (LD) analysis was carried out using HaploView 4.2 [49]. The Gene Structure Display Server 2.0 (https://gsds.gao-lab.org/) was used to visualize gene structure and SNP position. All plant materials were obtained from the seed bank of Yangtze University, College of Agriculture.

Statistical analysis

All experiments were conducted in three biological replicates. Data are expressed as means with standard deviation. Two-sided Student’s t-tests were performed to determine the statistical significant differences between means with p-values < 0.05 as threshold. The statistical analyses were conducted using Rstudio software (https://posit.co/products/open-source/rstudio/) (R version 3.6.1).

Results

Features of QTL for PHS/ seed dormancy in rice

Several QTL studies have been conducted to identify QTL regions associated with PHS resistance in rice, using diverse phenotypic characters/ traits to evaluate variations. Pre-harvest sprouting, seed dormancy, and low-temperature germination are the most frequently used traits (Table 1). Also, PHS resistance has been linked to other measurable agronomic traits such as GR, germination index (GI), and mean germination time (MGT). Herein, we compiled 167 PHS-associated traits QTL from 32 distinct mapping studies (Additional file 1: Table S2 and Table 1). The QTL were identified from second-filial generation populations (F₂), double-haploid, recombinant inbred lines, and backcross populations with population sizes ranging from 28 to 2,592. The QTL studies utilized various molecular markers, including SNPs, Simple Sequence Repeats (SSRs), Restriction Fragment Length Polymorphisms (RFLPs), Sequence Tag Sites (STSs), Cleaved Amplified Polymorphic Sequences (CAPS), Kompetitive Allele-specific PCR (KASP) and Functional Nucleotide Polymorphisms (FNPs) with markers ranging from 6 to 1,024 (Table 1). The number of QTL found in individual studies varied from 7 to 24, with chromosome 3 harboring the highest QTL (24), while chromosome 12 contained the lowest number of QTL (7) (Fig. 2a, b). The phenotypic variance explained (PVE or R²) by QTL varied from 0.06 to 89%, with most of the QTL possessing a PVE range of 6 to 15% (Fig. 2c). The LOD values of individual QTL ranged from 1.44 to 35.8 (Fig. 2d).

Fig. 2

Basic features of the initial QTL. a Distribution of QTL on 12 rice chromosomes, b number of QTL underlying different traits, c distribution of PVE values, and d, distribution of LOD scores

Identification of mQTL regions for PHS resistance in rice

We developed a consensus linkage map for all 12 rice chromosomes for QTL projection. The reference map consisted of 17,606 markers with a total genetic length of 2,299.05 cM (Additional file 1: Table S1). Out of the total initial QTL (167) collected, 134 (80.24%) were successfully projected onto the reference map (Additional file 1: Table S2). Markers on each chromosome were sparsely spread with opposite ends of the chromosomes displaying unequal marker density (Additional file 2: Fig. S1). Meta-analysis clustered these 134 QTL to 20 mQTL using Veyrieras algorithm in BioMercator v4.2.3 (Additional file 1: Table S3), of which three mQTL were singletons (Additional file 1: Table S3). The best model for individual chromosomes was chosen based on the values of AIC, AIC correction (AICc), AIC 3 candidate (AIC3), BIC and the AWE, respectively. The most projected QTL were found on chromosomes 7, 11, and 3, with chromosome 12 containing the lowest number of QTL projected (Fig. 3a). Moreover, most of the mQTL identified constitute initial QTL from diverse parental backgrounds and were from different environments suggesting their reliability and robustness in improving PHS resistance. These mQTL were unevenly distributed on all 12 chromosomes, ranging from 1 mQTL on chromosomes 5, 6, 9 and 10 to as many as 3 on chromosomes 7 and 11 (Additional file 1: Table S3; Fig. 4). The detected mQTL had mean LOD scores ranging from 2.65 to 14.08 and PVE² values ranged from 0.07 to 14.55% (Additional file 1: Table S3). The number of initial QTL involved in the formation of an mQTL varied from a singleton to 12. Of the 20 mQTL, 10 mQTL had at least 5 of the initial QTL (Additional file 1: Table S3). For CI, mQTL displayed a range of 0.24 cM to 7.61 cM distance, with a mean CI of 2.54 cM (Fig. 3b). Conversely, the projected QTL showed a wide CI variation ranging from 13.45 cM to 29.97 cM distance and a mean CI of 24.14 cM. As expected, the mQTL displayed a significant fold reduction in CI (9.56) compared to the initial projected QTL (Fig. 3b).

Fig. 3

Features of mQTL. a Histogram showing the distribution of initial QTL and the QTL projected, and b fold reduction in confidence interval (CI) of mQTL after meta-analysis

Fig. 4

Distribution of mQTL on12 rice chromosomes. mQTL regions are shown with the colored blocks inside each chromosome. Different colors on the left side of the maps indicate the initial QTL involved in mQTL. Names of markers and their positions are given on the right side of the chromosomes

Breeders’ mQTL for pre-harvest sprouting/ seed dormancy improvement

We implemented stringent criteria to select the most reliable and robust mQTL that could have potential implications in breeding programs for PHS resistance. These mQTL were selected based on the following criteria: (i) CI ≤ 2 cM, (ii) PVE ≥ 10, (iii) LOD ≥ 3, (iv) involvement of at least 5 initial QTL from diverse studies (Table 2). As a result, six of these mQTL were selected as “breeders’ mQTL” and recommended for use in breeding for PHS resistance in rice and subsequent analysis including candidate gene mining. The PVE score of breeders’ mQTL ranged from 10.42 to 14.55%, with the highest PVE observed for mQTL8.1, suggesting a significant contribution of breeders’ mQTL to phenotype variation explained (Table 2).

Table 2 Description of breeders’ mQTL associated with PHS/seed dormancy

Candidate gene mining within breeders’ mQTL regions and their functional attributes

The genes within the six breeders’ mQTL loci were considered promising CGs. Therefore, we performed batch retrieval of CGs from Rice Genome Annotation Project Database, resulting in 559 genes residing within the six mQTL intervals (Additional file 1: Table S4; Table 2). The genes were unevenly distributed in the breeders’ mQTL regions with the lowest (01) number of genes attributed to mQTL1.2; while mQTL10.1 contained the highest (408) number of genes (Table 2). In these mQTL, clusters of genes related to specific gene families were repeatedly identified which encode proteins belonging to different families including, zinc finger domain-containing protein, aminotransferase, leucine-rich repeats domain of unknown function (DUF), protein kinase domain, cytochrome P450 superfamily, no apical meristem (NAC), glutathione S-transferase, cellulose synthase-like protein, and GRAS family (Fig. 5a). To pinpoint potential CGs regulating PHS precisely, we incorporated DEGs transcriptome dataset for PHS [44] into CGs found in the breeders’ mQTL regions. Consequently, 34 common genes were identified between the two studies and thus were considered potential CGs for further analysis (Additional file 1: Table S5; Fig. 5b).

Fig. 5

Analysis of transcriptomic datasets and genes in breeders’ mQTL interval. a Frequency of candidate genes encoding known protein families associated with PHS, b Venn diagram of differentially expressed genes (DEGs) between transcriptome data and genes in the breeders’ mQTL interval, c GO terms for common genes in breeders’ mQTL regions, d KEGG enrichment pathways for common genes in breeders’ mQTL regions

Gene ontology analysis revealed the involvement of common CGs in diverse functional roles modulating PHS resistance including genes involved in cation binding, hydrolase activity, carbohydrate metabolic process, cellular protein modification process, and those associated with cell (Fig. 5c). Notably, key signaling pathways were highlighted for their possible functional relevance in regulating PHS resistance. The common CGs were involved in pathways related to carbon fixation in photosynthetic organisms, amino sugar and nucleotide sugar metabolism, arginine biosynthesis, and carbon metabolism (Fig. 5d).

In Silico expression of common CGs and protein-protein interaction of priori CGs

To identify potential CGs influencing PHS resistance in rice, we explored the expression pattern of the 34 common CGs using RNA-seq data from a public database (Rice Genome Annotation Project Database). The tissues included endosperm (25 DAP), seed, and embryo (25 DAP). The results revealed that LOC_Os10g18364, LOC_Os10g21930, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140 exhibited the highest expression levels in endosperm, seed, and embryo (Fig. 6). Consistent with the expression result, the eFP atlas showed high expression of these genes in various rice tissues including seed and inflorescence (Fig. 7c, d), which were closely associated with PHS. These results indicate that the CGs may play a vital role in PHS resistance in rice. Next, we investigated the potential interacting proteins of CGs using their protein sequences as query (Additional file 2: Fig. S2). Notably, LOC_Os10g18364 interacted with known gene ZEP1 which encodes synaptonemal complex protein 2 involved in chromosome synapsis and crossover during meiosis [80] and API5, an apoptosis inhibitor 5-like protein involved in the regulation of tapetal programmed cell death and degeneration during anther development [81]. Likewise, LOC_Os10g22590 a pollen Ole e I allergen and extension family protein interacted with three of its protein family members (LOC_Os05g27970, LOC_Os02g01190, and LOC_Os02g21280) and notably LOC_Os11g39370, a BRASSINOSTEROID INSENSITIVE 1-associated receptor kinase 1 precursor, LOC_Os11g41140, a zinc finger, C3HC4 type domain-containing protein, LOC_Os03g29410, a tyrosine-protein kinase domain-containing protein, and LOC_Os11g38020, GTPase of unknown function domain-containing protein (Additional file 2: Fig. S2). Moreover, LOC_Os10g25140, an aminotransferase domain-containing protein interacted with its family members with the same conserved domains including the Aminotran domain as LOC_Os05g39770, LOC_Os03g21960, LOC_Os06g35540, and LOC_Os03g07570 (Additional file 2: Fig. S2).

Fig. 6

Representative heatmap showing expression patterns of common genes in different rice tissues. All transcriptomic data were downloaded from Rice Genome Annotation Database. The change in colour signifies alteration in level of expression from low to high

Fig. 7

Relative expression profile of candidate genes and seed germination assay. a The seed dormancy phenotypes of 23xS-261 (resistant) and 23xS-262 (susceptible) seeds at 45 after heading (DAH). The phenotypes before and after germination assay, are presented. Scale bar, 10 mm, b Seed germination rates of 23xS-261 and 23xS-262 at 45 DAH. The significance was determined using the Student’s t-test, P < 0.05, c-d eFP rice atlas expression profile of selected candidate genes in different tissues of rice, e qRT-PCR analysis for selected candidate genes in rice embryo after 12 h of imbibition. Values are the means of three technical replicates ± standard error

Peak markers exhibit distinct polymorphism for breeders’ mQTL loci

To validate loci associated with PHS resistance in rice, two peak molecular markers linked to mQTL5.1 and mQTL10.1 were selected. The flanked markers, RG207 (mQTL5.1) and RM496 (mQTL10.1) exhibited distinct polymorphic fragments among the 16 rice genotypes evaluated (Additional file 2: Fig S3; Fig. 8a, b). Distinct DNA band patterns for marker RG207 with amplicon sizes of 220 bp and 240 bp among the genotypes were recorded (Fig. 8a). Interestingly, a similar trend was observed for SSR primer RM496 producing same band patterns with amplicon sizes of 220 bp and 240 bp among genotypes examined (Fig. 8b). These results indicate the presence of breeders’ mQTL loci suggesting that these loci could be validated further and introgress into elite background to enhance pre-harvest sprouting resistance in rice.

Fig. 8

Validation of mQTL loci controlling PHS/seed dormancy based on meta-QTL analysis. a Marker RG207, peak marker of mQTL5.1, b Marker RM496, peak marker of mQTL10.1. M: DL Marker 5000 bp; Number 1–16 represents the rice genotypes/cultivars showing band sample polymorphism of markers

To ascertain whether the observed polymorphism underlying the loci validated in genotypes exhibits a correlation with PHS, we specifically selected two rice genotypes 23xS-261 and 23xS-262 with distinct polymorphic fragments. Germination test under laboratory conditions showed that 23xS-261 line remained dormant with 28%, while 23xS-262 exhibited susceptibility to PHS with GR of 65% under relatively high humidity and high temperature (Fig. 7a-b).

qRT-PCR based verification of newly identified genes and phenotypic variation in pre-harvest resistant trait in rice

To verify the expression of newly identified genes for PHS, we compared the relative expression level between PHS susceptible line (23xS-261) and PHS resistant line (23xS-262) (Fig. 7e). The relative expression levels of four genes, namely LOC_Os10g18364, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140 were observed to be higher in PHS resistant genotype (23xS-261). Interestingly, LOC_Os10g25140 gene encoding aminotransferase domain-containing protein showed the most significant fold change expression value (2.23) in embryo in 23x-261 compared to 23x-262 followed by gene LOC_Os10g22590 with relative expression fold change of 1.60 in 23x-261 compared to 23x-262. The gene LOC_Os10g21940 which encodes an E3 ubiquitin ligase involved in syntaxin degradation exhibited a relative expression fold change of 1.50 in 23x-261 compared to 23x-262, while LOC_Os10g18364 displayed the lowest relative expression of 1.18 value in 23x-261 compared to 23x-262 (Fig. 7e).

A panel comprising 50 diverse rice cultivars was selected from an international rice panel. Germination rate (GR) was measured and used for association analysis. In detail, a wide variation in GR was observed, ranging from 47.5 to 100% (Additional file 1: Table S7). The GR distribution was continuous and normally distributed, indicating an extensive variation in our experimental population (Additional file 2: Fig. S4a, b).

Genetic variation in the intragenic regions of CG contributes to PHS in rice

To confirm the relationship between the four candidate loci and PHS, we obtained SNPs in intragenic regions of these genes using a panel of 50 accessions from the 3 K SNP-seek project. LOC_Os10g18364, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140 have intragenic variation underlying PHS; therefore, we further explored for superior haplotype alleles. Nine SNPs in LOC_Os10g18364 including seven intronic SNPs and two non-synonymous SNPs (missense variants) displayed significant association with PHS and were in complete linkage disequilibrium (LD) (Fig. 9a, b). The 50 accessions were grouped into three major haplotypes viz., Hap1, Hap2, and Hap3. A student t-test showed significant differences between the haplotypes (P = 0.0053) with Hap1 having the lowest mean GR (63.5%) (Fig. 9c).

Fig. 9

Natural variation of LOC_Os10g18364 and LOC_Os10g21940 alleles in cultivated rice. a,d. Survey of single nucleotide polymorphisms in intragenic regions ofLOC_Os10g18364 and LOC_Os10g21940, b,e. Linkage disequilibrium plot for the significant SNPs of LOC_Os10g18364 and LOC_Os10g21940, c,f. Boxplot showing GR in 50 rice cultivars as a function of each LOC_Os10g18364 and LOC_Os10g21940 haplotypes

The haplotype analysis of LOC_Os10g21940 showed two major haplotypes with Hap1 and Hap 2 with 7 and 43 accessions, respectively (Fig. 9). Six SNPs categorized these haplotypes with four SNPs (SNP11310024, SNP11310780, SNP11311506, and SNP11316533) being non-synonymous and exhibited complete linkage disequilibrium (Fig. 9d, e). Significant differences in GR were observed among these two major haplotypes with Hap1 recording GR of 63.6%, while Hap2 had 81.9% (Fig. 9f).

Three significant SNPs located within 5’UTR regions (SNP11723452) and exon regions (SNP11723540 and SNP11726136 (missense variants)) of LOC_Os10g22590 were associated with PHS and exhibited complete linkage disequilibrium (Fig. 10a, b). These SNPs grouped the accessions into three haplotypes, i.e., Hap1, 2, and 3 with 23, 25, and 2 accessions, respectively (Fig. 10c). The rice accessions with Hap1 displayed an average GR of (73.2%), which was significantly lower (P = 0.0154) than those accessions in Hap2 (85.4%) (Fig. 10c).

Fig. 10

Natural variation of LOC_Os10g22590 and LOC_Os10g25140 alleles in cultivated rice. a,d. Survey of single nucleotide polymorphisms in intragenic regions ofLOC_Os10g22590 and LOC_Os10g25140, b,e. Linkage disequilibrium plot for the significant SNPs of LOC_Os10g22590 and LOC_Os10g25140, c,f. Boxplot showing GR in 50 rice cultivars as a function of each LOC_Os10g22590 and LOC_Os10g25140 haplotypes

Furthermore, six variant sites (missense mutations) in the coding region of LOC_Os10g25140 were strongly associated with PHS and displayed complete linkage disequilibrium (Fig. 10d, e). These SNPs grouped the 50 accessions into three haplotypes and the GR of both Hap1 (64.1%) and Hap3 (70.3%) were significantly lower compared to Hap2 (85.5-100%) (Fig. 10f).

Overall, Hap1 and Hap3 associated with the four CGs emerge as distinct haplotypes underlying PHS. Therefore, we infer candidate genes LOC_Os10g18364, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140 to be associated with PHS in rice.

Discussion

Pre-harvest sprouting, primarily driven by unpredictable climate and untimely rainfall during seed maturation, severely impacts cereal crop productivity [82]. Breeding PHS-resistant rice varieties is imperative to alleviate the impact on global cereal yield and seed quality. Moreover, the efforts in breeding for PHS resistance have yielded limited results due to the complexity of the trait [83]. The identification of QTL associated with PHS has helped shed light on the regulatory mechanisms underpinning PHS in diverse crops including rice [79]. However, incorporating these QTL in breeding programs to improve PHS resistance often poses numerous challenges to breeders due to constraints associated with the detected QTL in various studies [22]. With these challenges in mind, mQTL analysis has been proposed to overcome these limitations [26]. Meta-analysis involves the integration of multiple QTL data from different studies carried out in various environments and genetic backgrounds, enabling the detection of stable and robust mQTL with reduced CI [41]. The present study represents the first attempt to identify mQTL governing PHS resistance in rice.

In this study, we performed a meta-analysis on 167 reported QTL governing PHS from 32 independent QTL studies, all focused-on mining genes/markers associated with PHS resistance in rice (Fig. 1). We projected a total of 132 QTL out of 167 QTL onto the reference map (Fig. 3). The rationale for the non-projection of the remaining QTL (35) comprising 20.96% of the initial QTL could be ascribed to non-existence of common markers or low marker density [39, 84]. The projected QTL (132 initial QTL) was refined into 20 mQTL based on QTL clusters. We found that the collated initial QTL (167) were unevenly distributed across the 12 chromosomes of rice, with chromosomes 3, 7, and 11 having the highest number of initial QTL (Fig. 1). However, chromosome 12 contained the lowest number of QTL reported (Fig. 2). Intriguingly, a greater majority of the projected QTL (approximately 55.30%) exhibited phenotypic variance > 10%, indicating that specific genetic regions might have impacted PHS resistance due to the sensitivity of the trait to environment [79]. Also, allelic variation at specific locus might confer stronger dormancy, hence the high PVE of initial QTL [85]. It is noteworthy that the initial QTL refinement in this study resulted in a significant fold reduction of CIs of the mQTL by 9.56 via meta-analysis (Fig. 3). Among the 20 mQTL detected, we found 10 mQTL that contained at least five initial QTL from different environments and diverse genetic backgrounds, indicating the reliability and robustness of the identified mQTL in this study [23, 86].

Furthermore, we implemented stringent criteria to select robust Breeders’ mQTL based on their low CI (≤ 2), high number of initial QTL involved in mQTL (at least 5), and high PVE values (> 10) (Table 2). Consequently, six Breeders’ mQTL namely mQTL1.2, mQTL2.1, mQTL5.1, mQTL8.1, mQTL10.1, and mQTL12.1 were identified (Table 2). All six Breeder’s mQTL were associated with multiple component traits suggesting the presence of pleiotropic genes or close interactions between genes that underpin different parameters [48]. The average CI of respective Breeders’ mQTL ranged from 0.07 to 1.98, and the average phenotypic variance explained ranged from 10.42 to 14.55%. Notably, mQTL8.1 had the highest PVE value among the Breeders’ mQTL (Table 2).

mQTL constitutes promising genomic loci that potentially contain high-confidence CGs underlying key agronomic traits. In this study, we extracted CGs within the six breeders’ mQTL, which resulted in 559 genes. We identified several genes to be repeated in various Breeders’ mQTL loci and on different chromosomes. The CGs found in multiple copies included proteins encoding zinc finger domain-containing protein, leucine-rich repeat, protein kinase, aminotransferase protein, no apical meristem protein, F-box domain-containing protein, amino acid transporters, cellulose synthase-like protein, GRAS family, cytochrome P450 superfamily, amino acid transporter, calcium-binding protein, and the most repeated copies being protein of unknown function (DUF) (Fig. 5a). In the previous study, seed dormancy 4 (SDR4) encodes a protein of unknown function (DUF) has been reported to be involved in the establishment of seed dormancy [10, 13]. Also, a weak seed dormancy 1 which encodes an aminotransferase protein revealed its crucial role in regulating seed dormancy and PHS in rice [87]. Notably, a cloned gene OsJAZ13 which encodes jasmonate ZIM-domain protein was located within Breeders’ mQTL10.1 region and plays vital roles in salt tolerance and modulating plant growth and cell death in rice [88, 89]. This adds credence to the consensus regions identified in this study. It is noteworthy that the presence of many DUF proteins within Breeders’ mQTL regions may be attributed to its potential involvement in signal transduction particularly in the ABA signaling pathway associated with PHS resistance in rice (Fig. 5a).

We deployed a robust and efficient strategy to pinpoint CGs underlying PHS resistance in rice by comparing genes in Breeders’ mQTL regions with publicly available differentially expressed transcriptome datasets. As a result, 34 common genes were found to overlap (Fig. 5b). Furthermore, we annotated the common genes (34) using GO and KEGG pathway analysis (Fig. 5candd) to gain insight into key biological processes and pathways involved in PHS resistance in rice. The common CGs were significantly enriched in carbohydrate metabolism process, cellular protein modification, cation binding, hydrolase activity, and cell. Surprisingly, we found that the common CGs were significantly enriched in key pathways including amino sugar and nucleotide sugar metabolism, carbon metabolism, and carbon fixation in photosynthetic organs, which have been implicated in PHS resistance in rice [90]. Next, we explored the expression profile of the common CGs using publicly available expression data in seed, embryo, and endosperm. Accordingly, we found five genes (LOC_Os10g18364, LOC_Os10g21930, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140) highly expressed in seed and embryo (Fig. 6). Consistently, developmental atlas expression data from the Rice eFP Browser showed that the five CGs were specifically highly expressed in seed and inflorescence which were closely related to PHS (Fig. 7). Therefore, we considered these genes potential candidate genes for further analysis. However, none of these genes have been previously characterized and functionally verified for their role in PHS resistance in rice. It is therefore possible to characterize these genes and validate their functions using a comparative genomics approach, followed by the development of functional markers.

We further subjected the candidate genes to expression analysis via qRT-PCR based on prior expression profile and their likely involvement in PHS resistance pathway in rice. The expression analysis was performed using two rice genotypes: 23x-261 (PHS susceptible line) and 23xS-262 (PHS resistant line) that exhibited significant differences in germination rate, and showed clear seed dormancy variation when subjected to high temperature and relatively high humidity conditions (Fig. 7). The expression result demonstrated increased fold expression of four priori genes (LOC_Os10g18364, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140) in PHS resistance genotypes than PHS susceptible genotype (Fig. 7). LOC_Os10g25140 encodes aminotransferase protein emerged as the gene with highest fold increase (2.23) in expression in PHS resistance material (23x-261) than in PHS susceptible genotype (23x-262). However, gene LOC_Os10g21930 which encodes F-box domain and kelch repeat-containing protein exhibited an increased fold expression in PHS susceptible genotype (23x-262) compared to PHS resistance line (23x-261) (Fig. 7). In previous studies, the overexpression of OsFbx352 enhanced GR of seeds when treated with ABA and decreases ABA biosynthesis genes (OsNced2 and OsNced3) indicating that F-box protein is mainly associated with seed germination rather than inhibiting germination to improve PHS resistance in rice [91, 92]. Seed dormancy and germination are influenced by internal and external cues [93]. The hormone ABA is largely known to function in inhibiting seed germination [94]. In line with this, most dormancy genes are involved in ABA synthesis and ABA signal transduction [13, 95, 96]. In the case of wheat, RsRHA2b, a gene that encodes E3 ubiquitin ligase, improves dormancy and PHS tolerance in transgenic wheat via ABA signaling pathway [97]. Moreover, MEDIATOR OF ABA-REGULATED DORMANCY 1 (MARD1), a zinc-finger gene in Arabidopsis mediated ABA-regulated seed dormancy [98]. Recently, a novel zinc-finger protein, OsZFP15, was discovered to act as a positive regulator of ABA catabolism. The overexpression of OsZFP15 decreased ABA content mainly through ABA catabolism and accelerated seed germination [99]. The pollen Ole e I allergen and extension family proteins have been recognized for their significant roles in pollen growth, cell wall metabolic activities, and stress response [100, 101]. However, their direct involvement in the ABA pathway, especially seed dormancy or PHS, is not well established. Nevertheless, their roles in stress response can be influenced by plant hormone signaling pathways such as ABA which require further study. These findings substantiate the potential roles of the newly identified genes in regulating PHS and therefore call for cloning and functional validation toward manipulating PHS in rice.

Protein-protein interaction network is crucial in understanding complex regulatory mechanisms involved in PHS and revealing potential targets for breeding or genetic engineering to improve PHS resistance in rice. For instance, the rice viviparous1 (Osvp1) gene interacted with transcription factors Rc and OsC1 to increase ABA sensitivity and improve PHS resistance [95]. Likewise, qSDR3.1 functions as a mediator of the ABA signaling pathway that regulates seed dormancy by interacting with ABSCISIC ACID INSENSITIVE5 (ABI5) [102]. In the present study, our protein-protein interaction prediction network revealed key proteins interacting with CGs (Additional file 2: Fig. S2). Notably, LOC_Os10g18364 encodes C2H2 zinc finger protein interacted with ZEP1 which encodes synaptonemal complex protein 2 involved in chromosome synapsis and crossover during meiosis [80] and API5, an apoptosis inhibitor 5-like protein involved in tapetal programmed cell death and degeneration during anther development [81]. The peak markers of Breeders’ mQTL5.1 and mQTL10.1 were validated using rice genotypes to ascertain their association with PHS or seed dormancy in rice. We found distinct polymorphism among the genotypes, indicating the existence of mQTL regions (Fig. 8). Importantly, it may be feasible to integrate the Breeders’ mQTL information in genomic prediction-based trait improvement programs by treating mQTL information as additive effects when making predictions using genomic data [103]. These markers have potential implications in marker-assisted selection (MAS) for breeding PHS resistance in rice.

The advent of biotechnological breeding has made significant progress in advancing the seed industry toward the realization of “Breeding 4.0”. Many studies have focused on genes associated with agronomic traits and emphasized elite natural variation applicable to breeding [104, 105). Using GWAS, Tai et al. [105] detected phosphatidylinositol 4-kinase (PI4K) family protein, TaPI4K-2 A, which is associated with PHS. Two distinct haplotypes were detected within TaPI4K-2 A, with Hap1 displaying lower germination and higher expression of TaPI4K-2 A in response to ABA sensitivity compared to Hap2 with the potential of integrating Hap1 in future wheat breeding programs to improve PHS resistance. In the present study, we obtained SNPs in intragenic regions of four priori genes to investigate functional natural variations and their impact on PHS. We identified a total of 10 major haplotypes in all four genes, with Hap1 belonging to LOC_Os10g18364, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140 exhibiting significantly lower germination rates (Figs. 9 and 10). For LOC_Os10g25140, we uncovered two favourable haplotypes Hap1 and Hap2 with a significantly lower germination rate. Non-synonymous variants specifically missense mutation were found in the coding regions of the four priori genes changing amino acid from one form to the other. A non-synonymous SNPs (G/C and G/A) were identified in the third exon of LOC_Os10g18364 which caused significant differences in GR among haplotypes. A total of four missense variant sites (G/A, T/C, G/A, and G/A) were identified in the coding region of LOC_Os10g21940contributing to PHS variation among haplotypes (Fig. 9). A study involving 263 stress-responsive genes in rice found that non-synonymous variants/SNPs in genes including MYB transcription factor family, WRKY transcription factor, calcium-dependent protein kinase, and sucrose transporter could affect encoded protein function and potentially influence drought and disease resistance in rice [106]. Apart from a missense mutation in the coding regions of LOC_Os10g22590, we also discovered a single variant site (G/A) in the 5’UTR region exhibiting a significant association with PHS. A recent study also showed that a base transversion from G to C in the 5’UTR region altered the expression of Sdt97 and confers semi-dwarfism in rice [107]. However, whether the SNP mutation residing in coding regions and 5’UTR region of these genes are the ones that cause variation in gene function warrants further study. Nevertheless, the identification of these loci provides new insights into the genetic basis of PHS and could serve as functional markers for breeding PHS resistance varieties in rice. Based on the result of this study, we proposed a plausible model to decipher the function of the newly identified genes in response to PHS (Fig. 11). Natural variation under the control of favorable haplotypes of LOC_Os10g25140^Hap1/3, LOC_Os10g18364^Hap1, LOC_Os10g21940^Hap1, and LOC_Os10g22590^Hap1 could regulate PHS via ABA biosynthesis and signaling pathway and unknown pathway in response to high temperature and relative high humidity conditions. However, this model remains to be tested in real time to ascertain the role of these genes in enhancing seed dormancy in rice.

Fig. 11

A proposed model of newly identified genes in PHS regulation in rice

Conclusions

The results of our meta-analysis yielded 20 mQTL associated with PHS of which six were designated as breeders’ mQTL based on narrow CI, high PVE, and number of initial QTL clusters. Candidate gene mining within the breeders’ mQTL loci led to the identification of 559 genes. Combined differentially expressed CGs from the public transcriptome dataset and 559 gene models revealed 34 common genes underlying PHS in rice. Furthermore, in silico expression analysis of the common genes pinpointed five genes with high expression in embryo and endosperm. The expression pattern of these CGs was further validated in PHS resistance and susceptible genotypes reaffirming their possible involvement in modulating PHS in rice. The natural variation identified in LOC_Os10g18364, LOC_Os10g21940, LOC_Os10g22590, and LOC_Os10g25140, particularly LOC_Os10g18364^Hap1, LOC_Os10g21940^Hap1, LOC_Os10g22590^Hap1, and LOC_Os10g25140^Hap1/Hap3 alleles hold promise for enhancing PHS resistance and represent a potential target for development of PHS resistance cultivars in rice.