Genome wide screening and comparative genome analysis for Meta-QTLs, ortho-MQTLs and candidate genes controlling yield and yield-related traits in rice

Background Improving yield and yield-related traits is the crucial goal in breeding programmes of cereals. Meta-QTL (MQTL) analysis discovers the most stable QTLs regardless of populations genetic background and field trial conditions and effectively narrows down the confidence interval (CI) for identification of candidate genes (CG) and markers development. Results A comprehensive MQTL analysis was implemented on 1052 QTLs reported for yield (YLD), grain weight (GW), heading date (HD), plant height (PH) and tiller number (TN) in 122 rice populations evaluated under normal condition from 1996 to 2019. Consequently, these QTLs were confined into 114 MQTLs and the average CI was reduced up to 3.5 folds in compare to the mean CI of the original QTLs with an average of 4.85 cM CI in the resulted MQTLs. Among them, 27 MQTLs with at least five initial QTLs from independent studies were considered as the most stable QTLs over different field trials and genetic backgrounds. Furthermore, several known and novel CGs were detected in the high confident MQTLs intervals. The genomic distribution of MQTLs indicated the highest density at subtelomeric chromosomal regions. Using the advantage of synteny and comparative genomics analysis, 11 and 15 ortho-MQTLs were identified at co-linear regions between rice with barley and maize, respectively. In addition, comparing resulted MQTLs with GWAS studies led to identification of eighteen common significant chromosomal regions controlling the evaluated traits. Conclusion This comprehensive analysis defines a genome wide landscape on the most stable loci associated with reliable genetic markers and CGs for yield and yield-related traits in rice. Our findings showed that some of these information are transferable to other cereals that lead to improvement of their breeding programs.

as quantitative trait loci (QTLs) [2,9], dealing with them is a challenge. QTL mapping provides accurate deciphering of genomic regions regulating these complex traits [10] and it has accelerated the success of breeders for improving quantitative traits by marker-assisted selection (MAS) [11]. However, the main problem faced by researchers in using QTL results are their dependency upon the population genetic backgrounds and the phenotyping environment that limit their applications in a wider range of populations or environments [10,12].
Meta-analysis of QTLs unravels consensus and stable QTLs by merging different QTLs from independent experiments regardless of their genetic backgrounds, population types, evaluated locations and years [12][13][14]. Therefore, the Meta-QTL, with the abbreviation of "MQTL" in the rest of the manuscript, results are highly reliable and they can be widely used in breeding programs. Moreover, MQTL analysis consistently refines the position of QTLs and narrows down the confidence intervals (CI) that leads to accuracy of MAS [15,16]. This conceptual approach has been used to detect MQTLs for various traits in barley [16,17], wheat [11,[18][19][20], soybean [21,22] and maize [15,[23][24][25][26][27]. In rice, there are two MQTL studies on YLD, PH and TN traits. Of these, one was conducted on 11 QTL studies published from 1998 to 2008 [28], whereas another was performed on 35 QTL studies covering the period of 1995 to 2006 [29]. Moreover, Daware et al. (2017) reported seven MQTLs related to GW from 7 QTL studies published only since 2008 to 2015 on indica and aromatic rice accessions [10].
We conducted a large and comprehensive metaanalysis on QTLs of YLD, TN, GW, PH and HD traits that are reported from 101 studies published from 1996 to 2018 in 122 bi-parental populations evaluated under unstressed conditions. It is the most comprehensive MQTL study for aforementioned traits in cereals and the first MQTL study on HD in rice. Beside MQTL study, each of the detected MQTLs was investigated to identify candidate genes (CGs) related to the evaluated traits. In addition, due to high synteny among rice, barley and maize [30,31], we expanded our analysis to detect ortho-MQTLs in among these cereals. The uncovered novel MQTLs, ortho-MQTLs and candidate genes will aid genetic dissection of yield-related traits to improve yield in cereals.

Main features of yield-related QTL studies in rice
A total of 1052 QTLs controlling YLD, GW, HD, PH and TN in rice under unstressed conditions were retrieved from 122 populations reported in 101 studies since 1996 (Table 1). The number of QTLs for each trait and their distribution on 12 chromosomes of rice are presented in Fig. 1a and b. The QTLs scattered unevenly on different chromosomes; while chromosome 3 harbored the largest number of QTLs with 180 QTLs, followed by chromosome 1 (153 QTLs) and 7 (111 QTLs), chromosome 9 had the lowest number of QTLs with 36 QTLs.
The number of QTLs was varied in different evaluated quantitative traits. Among the studied traits, GW and HD had the highest number of QTLs with 339 and 267 QTLs, respectively, followed by PH, YLD and TN with 204, 165 and 77 QTLs, respectively (Fig. 1b). The QTLs for GW were mainly located on chromosome 3, 5 and 1 with 60, 48 and 48 QTLs, respectively, and the majority of QTLs for HD were placed on chromosomes 3 (56), 7 (44) and 6 (43). Consistently with previous reports [28,130], chromosome 1 had the highest number of QTLs for YLD. Chromosome 1 also harbored the highest number of QTLs for PH and TN traits (Fig. 1b).

Detected MQTLs for yield-related traits
A total of 960 QTLs out of the 1052 QTLs (91%) from 122 populations were successfully projected on the reference map ( Table 2). The MQTL analysis confined these QTLs into 114 MQTLs (11.87 %) with QTLs originated from at least two studies for all the aforementioned traits (Table 3; Fig. 1, 2 and S1). Of these MQTLs, 58 MQTLs (50.8 %) were obtained from at least three independent studies (Table 3; Additional file 1).
The number of MQTLs for each trait was distributed unevenly among rice chromosomes. In this analysis 34, 23, 28, 19 and 10 MQTLs were detected for GW, HD, PH, YLD and TN traits, respectively. The distribution of MQTLs for each trait on each chromosome is presented in table 3 and Additional file 1. The most of the MQTLs associated with GW were located on chromosomes 1 and 5, whereas MQTLs of HD were mainly located on chromosomes 3 and 7 (Table 3). Overall, we could detect at least one MQTL for GW on all of the chromosomes (Table 3). Apparently, chromosome 1 was predominantly involved in controlling PH, YLD and TN traits. The lowest MQTLs for GW, HD, PH, YLD and TN were mainly located on chromosomes 5, 9, 10, 11 and 12. In general, there was a positive correlation between QTLs density and the number of MQTLs on chromosomes for all studied traits (r=0.90, Table 2 and 3, Fig. 1b). Moreover, the traits with the higher number of QTLs had the higher number of MQTLs (Fig. 1a).
A MQTL with the higher number of initial QTLs is a more stable MQTL independent from genetic background and environment. MQTL-HD8 with 13 initial QTLs had the highest number of QTLs derived from 11 different populations followed by MQTL-HD5, MQTL-GW6 and MQTL-GW16 with 11, 10 and 10 initial QTLs derived from 11, 7 and 4 different populations, respectively (Table 3). These MQTLs appeared as the most robust, viable and stable QTLs in different locations and    years. Furthermore, we identified 22 overlapping MQTLs or clusters of MQTLs which controlled at least two traits (Additional file 1). Interestingly, two clusters of MQTLs located on chromosomes 7 and 8 includes all studied traits (Additional file 1). The overlapping MQTLs are likely to contain CGs with broad pleiotropic effects.
The distribution pattern of MQTLs on the rice genome was investigated and compared with genomic events including selective sweep regions and gene density. The number of MQTLs per chromosome varied from 2 (chromosome 12) to 21 (chromosome 1) with an average of 9.5 MQTLs per chromosome (Table 3; Fig. 2 and additional file 1). The overview on the distribution of gene density on the rice genome revealed that sub-telomeric regions harbor most of the genes ( Fig. 2 and 3). Similarly, the distribution of QTLs and MQTLs displayed comparable pattern to the gene density over the rice genome ( Fig. 2 and 3). We detected the lowest number QTLs at the centromeric intervals for all studied traits ( Fig. 2 and 3).
A total of 23 and 12 MQTLs were co-located on the selective sweep regions and the regions containing known functional variants on the rice genome, respectively [131]. These regions can be further investigated among the rice genetic resources for improving yield in breeding programs (Fig. 2, Additional file 3).
Beside known genes, we identified novel CGs based on their annotated function that are presented in Additional file 2 and potentially can be a regulator of GW. In MQTL-GW6 on    chromosome 2, the Os02g0283800 annotated as OsBAK1-5 or OsSERK2 regulates grain size and number [136]. The OsALMT7 gene located at MQTL-GW7 interval was shown to affect grain size [137]. In MQTL-GW10 interval on chromosome 3, OsEZ1 and NRL2 genes contribute to the grain size and YLD in rice [138,139]. Furthermore, MQTL-GW20 with six GW QTLs (Table 3) from three populations, contains the Os05g0414700 gene encoding a BRASSINOS-TEROID INSENSITIVE 1-associated receptor kinase 1 that might be a new CG for regulating GW. We also detected new gene which encodes brassinosteroid receptor kinase gene at MQTL-GW28 intervals that was recently shown to collocated at selective sweep regions and have high impact on development and therefore could be good candidates for further functional investigation.

MQTLs and CGs for Heading Date
HD controlled by polygenes that substantially encompasses the seed production and YLD [7]. In this analysis, we could detect 23 MQTLs related to HD including four MQTLs on chromosomes 3 and 7, three MQTLs on chromosomes 2 and 8, two MQTLs on chromosomes 4, 5 and 6 and one MQTL on chromosomes 1, 10 and 11.

MQTLs and CGs for Yield
YLD is the most prominent criteria in the rice breeding programs [96]. We identified 19 MQTLs for YLD comprising five MQTLs on chromosome 1, four MQTLs on chromosome 3, two MQTLs on chromosomes 2, 4 and 8, one MQTL on chromosomes 6, 7, 9 and 11 ( Table 3). The proved genes such as GIF2 [159], OsLSK1 [160], d11 [154], APO1 [161]/DEP3 [9] and OsSPL13 [135] controlling yield were located on MQTL-YLD3, YLD4, YLD13 YLD14 and YLD15 respectively, and most of them were co-located with selective sweep regions and functional variants on coding regions (Additional file 3) . The list of CGs is presented in Additional file 2. For instance, Os01g0171000 gene at MQTL-YLD1 which encodes BRASSINOSTEROID IN-SENSITIVE 1-associated receptor kinase is a potential candidate for higher yield. In MQTL-YLD6, the SWEET15 gene was reported to have substantially effects on YLD through regulating seed filling in rice [162].

MQTLs and CGs for number of tillers
TN is a foremost feature in plant architecture and grain production in cereals. Despite its agronomic importance, only a few tiller controlling genes have been identified so far [163,164]. A total of 10 MQTLs associated with TN are detected in our analysis including two MQTLs on chromosomes 1 and 3 and one MQTL on chromosomes 2, 4, 6, 7, 8 and 9 (  OsPIN5b gene at MQTL-TN10 on chromosome 9 controls TN, PH, panicle size and other traits related to plant architecture [167]. Moreover, the GSK3/SHAGGYlike kinase, GA2ox genes on MQTL-TN1 and GRAS TFs on MQTL-TN4 and MQTL-TN7 are also reported to regulate TN [163].

Ortho-MQTL mining in barley and maize
Due to the high synteny among barley, maize and rice and also the economically importance of the studied traits in all cereals, 58 of the most reliable MQTLs derived from at least three independent studies in rice were selected for investigation of ortho-MQTLs in barley and maize. Consequently, a total of 11 ortho-MQTLs were detected for rice and barley including four ortho-MQTLs for HD, three for GW and PH and one for YLD. Moreover, a total of 15 ortho-MQTLs were identified in rice and maize consisting of nine and six ortho-MQTLs for YLD and PH, respectively. Among them, three ortho-MQTLs (ortho-MQTL-PH6, ortho-MQTL-PH10 and ortho-MQTL-YLD11) were cross-species in all the three crops (Table 4; Fig. 4).
Moreover, three rice MQTLs (MQTL-PH6, MQTL-PH10 and MQTL-YLD11) on chromosomes 1, 2 and 3, respectively, were located in a syntenic region with detected MQTLs on barley chromosomes 3H, 6H and 4H, respectively, as well as maize MQTLs for YLD (Fig. 4a, b and c; Table 4). Ortho-MQTL mining could validate our analysis and it can facilitate detecting underlying regulatory genes with evolutionary history and conservative function. All the genes located at the ortho-MQTLs regions along with their annotations were reported in Additional file 4. Remarkably, the orthologous of wellknown proved genes including Hd6, OsbZIP62 and OsSPL13 in rice were detected in barley and maize ortho-MQTLs.

Distribution pattern of QTLs and MQTLs and identification of CGs
The investigated QTLs were not evenly distributed on all chromosomes of rice. Chromosomes 1, 3 and 7 harbor the largest number of QTLs that are in agreement with previous reports [2,168]. Chromosomes 3, 1 and 5 harbored the largest number of QTLs associated with GW in parallel with preliminary results in rice [2]. Additionally, HD had the highest initial QTLs on chromosomes 3, 7 and 6 as shown by prior results [7].
The MQTL analysis detects the most stable QTLs regardless of the genetic background, the phenotyping variations in different places and years, and the markers density that are the main restrictions of QTL mapping [8,[10][11][12]17]. Association mapping is another approach with higher accuracy in compare to QTL mapping for identification of genomic regions underlying quantitative traits [10]. But, this approach faces remarkable falsepositive results due to the structure of population used in the analysis [10]. Therefore, MQTL analysis is considered as the most reliable approach to identify stable loci controlling quantitative traits. Here, MQTL analysis confined a total of 960 QTLs into 114 MQTLs on twelve chromosomes of rice for all studied traits. Chromosomes 1 and 12 with 21 and 2 MQTLs had the highest and the least number of MQTLs, respectively. Similarly, Swamy and Sarla (2011) reported the highest number of MQTLs for YLD on chromosomes 1, 2 and 3 [28]. Beside the common MQTLs reported in other studies [10,28,29] for YLD, PH, TN and GW (Table 3), we identified 31 new MQTLs for GW, 21 new MQTLs for PH, 9 new MQTLs for TN and 7 new MQTLs for YLD and this is the first MQTL study conducted on HD (Table 3).
It is hypothesized QTL density is chiefly related to gene density and polymorphism rate [15]. Our results demonstrated that the most of MQTLs and QTLs were located at the subtelomeric regions where gene density is relatively high. Previous investigations in barley and maize reported similar results in which QTLs and MQTLs were densely located at the subtelomeric regions [15,17,169].
The higher sequence polymorphism rate and functional variants at coding regions resulted in higher differentiation in allele frequencies in rice [131]. We detected 12 MQTLs which were precisely located at these regions. In addition, 23 MQTLs were cope with the selective sweep regions that occurred during the domestication processes and profoundly affected the selection and spread of critical traits [131]. These results provide beneficial information for breeders to proficiently apply diverse genetic resources for improving rice and other cereal crops. We detected fundamental genes for GW, HD PH, and YLD (Additional file 3) controlling aforementioned traits in these overlapping MQTLs.
Moreover, MQTL analysis considerably reduces the CI in compared to the initial QTLs. It consequently diminishes the number of genes anchoring at the QTL interval that are in agreement with previous reports [16,23,27,130,168]. Therefore, MQTL analysis enhances the precision of CGs prediction and the detection of markers for  [24] marker assistant selection in breeding [17]. In our analysis, the average CI was reduced up to 3.5 folds in compare to the mean CI of the original QTLs, therefore, the number of genes located at the QTLs interval was extensively reduced. Fundamental genes such as Hd1, Hd5, Hd6, Hd17, DTH3, HDR1, OsMADS3, OsMDAS6, OsMADS18 and OsMADS22 for HD, and d2, Gn1a, d11, GS2, RSR1, GS5, OsSPL13 and SRS5 for GW were still located at the narrowed MQTLs interval. Furthermore, among the limited number of genes annotated at each MQTL interval we detected potential CGs for PH, YLD and TN attributes that are listed in Additional file 2.

Ortho-MQTL mining
Despite the high interest in identification of genes involved in YLD and yield-related traits in barley and maize as two economically important crops, the responsible genes have largely remained unknown due to their complex genomes. Given a close evolutionary relation among grass genomes [170], synteny analysis of barley and maize with rice as a model crop in grasses enabled us to broaden our genetic information among these species [30]. Identification of ortho-MQTLs among these close species expands their utility and it also validates their stability and the confidence of related CGs. Here we selected the most prospective rice MQTLs containing at least three QTLs from different studies to explore their conserved syntenic regions reported in similar MQTLs studies on the same traits in barley and maize to identify ortho-MQTLs (Table 4).
For rice MQTL-HD8 there is a MQTL in the syntenic region on barley [17] controlling ortho-MQTL-HD8 containing a rice Hd6 orthologous gene (HOR-VU5Hr1G097230) known to have a high impact on HD [7]. Moreover, in the syntenic region of rice MQTL-HD18 in barley there is a MQTL on chromosome 7 (ortho-MQTL-HD18) encompassing OsbZIP62 orthologous gene (HORVU2Hr1G017020) that regulates flowering in rice [143]. All other ortho-MQTLs identified between rice and barley are described in table 4 and the orthologous genes are presented in Additional file 4 among which of the most important CGs were explained during the corresponding rice MQTL discussion.

Comparison of MQTLs with GWAS studies
Detected rice MQTLs were compared with GWAS studies in rice and barley led to identification of common significant loci that provides more confident MQTLs. Consequently, in rice 7 and 11 significant GWAS signals for GW and HD, respectively, were co-located with our MQTLs (Additional file 5). They were distributed on all chromosomes of rice except chromosomes 9, 11 and 12. These results indicate the compatibility and coherence of the two methods in identification of significant genomic regions corresponding to the studied traits. Four rice MQTLs including MQTL-GW4, MQTL-GW20, MQTL-GW24 and MQTL-HD4 on chromosome 1, 5, 7 and 2, respectively, were in co-linear with the syntenic regions containing significant GWAS signals for the same traits in barley [172,173]. Five genes with a significant signal in barley GWAS had orthologous genes in the detected rice MQTLs intervals (Table 5), including Vrs3 gene on chromosome 1H [172] that its orthologous in rice located at MQTL-GW20 on chromosome 5, a  Table 4 and Table S3.
serine carboxypeptidase gene on barley 2H [173] and its orthologous in rice MQTL-GW24 on chromosome 7 that regulates GW [132]. However, this approach has limitations related to compatibility of the genome intervals between MQTLs and significant peaks in GWAS.

Conclusions
In conclusion, we succeeded to define a genome wide landscape on the most stable loci that associate with genetic markers and CGs related to yield and yieldrelated traits in rice. Our findings show that MQTLs of evaluated important agronomic criteria appear at least partly to be transferable to other cereals and genome wide association studies that helps breeding programs in cereals.  Table 1.

QTLs projection on reference map
The most comprehensive available genetic map integrated from six identified well-known genetic maps in rice, with the highest number and various types of markers [168], was selected as the reference map in the present study. This map is highly saturated with 6969 markers of 0.25 cM average distance that results in a total length of 1771.8 cM on 12 chromosomes with an average chromosome length of 147.65 cM. The position, chromosome groups, the proportion of phenotype variance (R 2 ), and the log of odds ratio (LOD score) were collected for each of the QTLs in the 122 used populations. In order to calculate 95% of CI for QTLs, we used the formula, CI=530/(N*R 2 ) for BC and F 2 lines, CI=287/(N*R 2 ) for DH lines and CI=163/(N*R 2 ) for RILs lines [174], where N is the population size and R 2 is the proportion of phenotypic variance of the QTL. For QTLs without precisely defined LOD scores [39,48,77] and R 2 [63,69,91], those criteria were arbitrarily quantified as 3 and 10%, respectively. All of the collected QTLs with proper information were projected onto the reference map by BioMercator V4.2 [12,14]. In order to be able to incorporate QTLs derived from the studies based on SNP markers in our MQTL analysis, the position of the flanking markers were ascertained on the rice genome and closest markers on the reference map were used in our analysis. Consequently, 960 out of 1052 initial QTLs were successfully projected on the reference map.

MQTL analysis
The meta-analysis was carried out for the integrated QTLs from different studies and relocated on consensus position of each MQTL using BioMercator V4.2 [12,14]. The algorithms and statistical procedures implemented in this software are well-described in the literature [12,14,175]. The best model of estimated MQTLs was selected based on the prevalent value among AIC (Akaike information content), AICc (AIC correction), AIC3 (AIC 3 candidate models), BIC (Bayesian information criterion) and AWE (average weight of evidence) The rice gene ID located at the corresponding MQTL interval with a barley orthologous gene containing a SNP with significant signal in GWAS analysis.
criteria and it was considered as the best fit. Consequently, the consensus QTL from the optimal model was reported as a MQTL. Mapchart V.2.32 software [176] was applied to demonstrate the MQTLs and related QTLs on the reference map. The position of MQTLs on the rice genome was shown as a heatmap using pheatmap and RIdeogram R package [177,178].
To investigate the distribution of MQTLs towards centromeric and telomeric regions in the genomic point of view, the centromere position was retrieved from Cheng et al. (2002) and Kawahara et al. (2013) studies [179,180] and shown on each chromosome. Moreover, to expand our genomic approaches, the distribution of MQTLs were compared to the selective sweep regions (domestication loci) and functional variants in coding regions with strong alteration in allele frequency between cultivated and wild rice [131]. Additionally, the distributions of gene density, QTLs and MQTLs were drawn on chromosomes using SOFIA R package [181]. Finally, all the detected rice MQTLs were compared with the position of significant loci related to the traits resulted from Genome Wide Association Studies (GWAS) using the Rice SNP-seek database [182].

Identification of candidate genes
To determine CGs related to YLD, TN, GW, PH and HD traits located at the corresponding region of each detected MQTL, the rice genome (IRGSP-1.0) was investigated in EnsemblPlants (https://plants.ensembl.org/ index.html) using the position of flanking markers obtained from the Gramene (http://archive.gramene.org/ qtl/) database. For those flanking markers without genomic position, the closest markers from consensus genetic reference map were exploited to project the MQTL on the genome. Consequently, all the genes underlying the genomic region of each MQTL were functionally annotated by EnsemblPlants and FunRiceGenes (https:// funricegenes.github.io/) [183], and CGs were introduced based on their description and putative function in rice and closely related species.