Pinpointing genomic regions associated with root system architecture in rice through an integrative meta-analysis approach

Applying an integrated meta-analysis approach led to identification of meta-QTLs/ candidate genes associated with rice root system architecture, which can be used in MQTL-assisted breeding/ genetic engineering of root traits. Root system architecture (RSA) is an important factor for facilitating water and nutrient uptake from deep soils and adaptation to drought stress conditions. In the present research, an integrated meta-analysis approach was employed to find candidate genes and genomic regions involved in rice RSA traits. A whole-genome meta-analysis was performed for 425 initial QTLs reported in 34 independent experiments controlling RSA traits under control and drought stress conditions in the previous twenty years. Sixty-four consensus meta-QTLs (MQTLs) were detected, unevenly distributed on twelve rice chromosomes. The confidence interval (CI) of the identified MQTLs was obtained as 0.11–14.23 cM with an average of 3.79 cM, which was 3.88 times narrower than the mean CI of the original QTLs. Interestingly, 52 MQTLs were co-located with SNP peak positions reported in rice genome-wide association studies (GWAS) for root morphological traits. The genes located in these RSA-related MQTLs were detected and explored to find the drought-responsive genes in the rice root based on the RNA-seq and microarray data. Multiple RSA and drought tolerance-associated genes were found in the MQTLs including the genes involved in auxin biosynthesis or signaling (e.g. YUCCA, WOX, AUX/IAA, ARF), root angle (DRO1-related genes), lateral root development (e.g. DSR, WRKY), root diameter (e.g. OsNAC5), plant cell wall (e.g. EXPA), and lignification (e.g. C4H, PAL, PRX and CAD). The genes located within both the SNP peak positions and the QTL-overview peaks for RSA are suggested as novel candidate genes for further functional analysis. The promising candidate genes and MQTLs can be used as basis for genetic engineering and MQTL-assisted breeding of root phenotypes to improve yield potential, stability and performance in a water-stressed environment.


Introduction
Rice (Oryza sativa L.) is one of the most important staple crops, which is consumed by more than one-third of the world's population. Among the abiotic factors, drought is the major limiting factor affecting rice growth and productivity. It is estimated that by the year 2050, more than 50% of the world's arable land will be affected by drought (Singhal et al. 2016). Therefore, improving rice for drought-prone environments is a priority.
Dynamic responses of root system architecture (RSA) play a key role in efficiently using water and nutrients in crops (De Smet et al. 2012;Paez-Garcia et al. 2015), and also in drought avoidance mechanism and drought tolerance performance (Blum 2005). Improving the architecture and function of roots will be the key to the second green Communicated by Matthias Wissuwa. revolution in the future (Gewin 2010). However, breeding for root morphological traits has rarely been addressed, mostly because of the costs and time constraints as well as lack of reliable and efficient phenotyping techniques for root-related traits (Carvalho et al. 2014;Toyofuku et al. 2015). The RSA indicates the organization of the primary roots, crown roots, and lateral roots in the soil environment (De Smet et al. 2012). The genetic factors and interactions with environmental conditions mainly affect the RSA adaptation (Malamy 2005). Genotypic background determines the intrinsic morphological characters, whereas environmental factors modify the root morphologies on the basis of requirements for dynamically adapting to undesirable environmental conditions (Gruber et al. 2013;Soriano and Alvaro 2019). Mapping of quantitative trait loci (QTL) has been used as a powerful statistical method to identify genomic regions associated with traits of interest for breeding (Wang et al. 2019b). To date, many QTLs related to RSA traits were identified by linkage analysis from different population types and sizes across diverse moisture conditions (Courtois et al. 2003;de Dorlodot et al. 2007;Uga et al. 2011Uga et al. , 2013. In many QTL mapping studies, the overlap of QTLs for root and agronomic traits has been revealed, suggesting the profound implications that breeding for RSA will have on improving crop genotypes through enhancing crop productivity and high water/nutrients use efficiency under water-limited conditions (Jia et al. 2019;Ju et al. 2018;Maccaferri et al. 2016;Tuberosa et al. 2002). The growth angle and length of roots mainly determine the root system distribution of rice throughout the soil profile (Araki et al. 2002).
Various positions of a QTL in different mapping populations lead to a large confidence interval and uncertain position of the QTL. This can be further exacerbated by several additional factors such as different sizes of mapping populations and sampling errors (Darvasi and Soller 1997;Darvasi et al. 1993), differences in experimental replicates, marker density and QTL mapping models. Various approaches have been used so far to confirm the QTL results, such as QTL mapping using first-generation populations, and getting validated in advanced generation breeding populations of the same cross (Gelli et al. 2017). In other cases, QTL confirmation is performed using the candidate gene approach or positional cloning, followed by integrating functional and genetic data within a breeding process (de Dorlodot et al. 2007). This challenging process, however, requires highdensity linkage maps, extensive genomic resources and logical informatics data (de Dorlodot et al. 2007).
To identify consensus QTL regions across multiple studies, QTL meta-analysis method was initially developed (Goffinet and Gerber 2000) using maximum likelihood estimation and was then improved (Veyrieras et al. 2007). Meta-analysis clearly estimates the numbers, positions and CI of meta-QTL regions in each chromosome. This method has been used to identify consensus regions of the genome across multiple QTL studies for their effect and consistency across different genetic backgrounds and environments, also to refine and confirm QTL positions on a consensus map via mathematical models. The QTL meta-analysis has been performed on root morphological traits in different species, including maize (Guo et al. 2018), bread wheat (Darzi-Ramandi et al. 2017;Soriano and Alvaro 2019;Bilgrami et al. 2020), durum wheat (Iannucci et al. 2017), and oilseed rape ).
In the current study, we integrated 425 QTLs from the 34 published QTL mapping studies for root traits in rice and identified consensus genomic regions through QTL metaanalysis method into a novel integrated consensus genetic map with 5447 loci. Thus far, only one QTL meta-analysis on rice root morphological traits has been reported based on the QTLs collected from 24 studies from 1995 to 2007, mostly using AFLP and RFLP markers, detecting MQTLs with 1.98 Mb average CI and a reduction of 61% in the number of QTLs (Courtois et al. 2009). However, out of 34 QTL mapping studies used in the current research, 27 studies, including 268 QTLs, were reported after 2007, using mostly SSR and some SNP markers, which had not been included in the QTL meta-analysis so far. Moreover, the supporting intervals of the identified MQTLs reduced to an average of 1.57 Mb corresponded to a reduction of 85% in the number of initial QTLs and then validated by the related GWAS studies. The genes located in the MQTLs were detected and functionally classified. In addition, the differentially expressed genes (DEGs) in the rice root under drought conditions were detected through the analysis of RNA-seq and microarray datasets, and the MQTL regions associated with RSA were explored to find the drought-responsive genes in the rice root. Conclusively, integrating QTLs, GWAS, and transcriptome data resulted in identification of the promising MQTLs and candidate genes, which can be used in MQTLassisted breeding and genetic engineering to improve RSArelated traits in rice.

Collecting RSA trait-associated QTLs from independent studies
A comprehensive review of the literature was conducted by searching the articles published from 2001 to 2020 on rice QTLs associated with RSA traits in normal conditions and drought stress ( Table 1). The data collected for each QTL region included (1) root architecture traits, (2) parents of the population, (3) types of mapping population: F 2 , backcross (BC), doubled haploids (DH), recombinant inbred lines (RILs) and near-isogenic lines (NILs), (4) population size (N), (5) logarithm of odds ratio (LOD score), (6) proportion of phenotypic variance explained by the QTL (R 2 ), (7) the flanking or single marker(s) for interval mapping and single-marker analysis (SMA), respectively. The present study assumed a LOD score of 3 in a few cases in which the published article reported a p-value statistic or stated that a minimum LOD score of 3 was adopted as the threshold for QTL analyses. Equation (1) proposed by (Nagelkerke 1991) was used to estimate the LOD score (where N represents the population size), which is the explanatory power of QTLs if the LOD score was not reported.
Ten traits associated with root architecture analyzed in the present study under normal and drought stress conditions included deep root ratio, root number (RN), root length (RL), root thickness (RT), root volume (RV), root to shoot ratio (RSR), root fresh weight (RFW), root dry weight 1 3 (RDW), deep root ratio (DRR), root growth rate (RGR) and root surface area (RSA) were analyzed. The analysis encompassed the QTLs whose map positions, LOD scores and R 2 values were available. All the QTL studies related to RSA conducted on these traits using markers including AFLP, SSR, SNP and RFLP were applied in the QTL meta-analysis, but the articles lacking proper genetic maps or QTL-associated data were excluded. The QTL mapping studies with any missing parameters were also discarded. All the markers from the aforementioned maps were used to develop the consensus marker map. The meta-analysis was performed after eliminating the markers showing inversion on the consensus map. The projection of the QTLs position was performed on the basis of a simple scaling method between the interval of the QTL flanking markers on their original map and the interval of these markers on the consensus map (Supplementary Table S1). After estimating the new confidence interval of the initial QTL on their original genetic map using the Gaussian distribution, the QTLs were projected on the consensus map. The 95% Cl for a QTL position was calculated according to the type of population used in Eqs. (2)-(4), in which N represents the population size, R 2 the percentage of the phenotypic variance of the QTLs and CI the supporting or confidence interval for each initial QTL.

Developing the consensus map and projecting the QTLs
( (2) is presented by modeling to calculate the 95% CI of F 2 and BC populations by (Darvasi and Soller 1997). Equation (3) was used for the population of DH lines (Visscher and Goddard 2004), and Eq. (4) for both RIL and NIL populations (Guo et al. 2006).

Meta-analysis and QTL-overview index
After projecting the QTLs on the consensus map, meta-analysis was conducted according to the QTL clusters on each chromosome using BioMercator V4.2 (Arcade et al. 2004) which contains algorithms from the MetaQTL software (Sosnowski et al. 2012;Veyrieras et al. 2007). Two different approaches were used depending on the number of initial QTLs in a chromosome, which was below 10 when the metaanalysis proposed by Gerber and Goffinet was employed. Based on this method, the most likely assumption is estimated in BioMercator among five MQTL models (1, 2, 3, 4, or N) with different AIC values. Moreover, the model with the minimum AIC was selected for QTL integration and identification of consensus MQTL positions. The method proposed by Veyrieras et al. was employed in case the number of QTLs in a chromosome equaled at least 10. According to this approach, two stages are taken in QTL metaanalyses. In the first stage, the collected QTLs are clustered on individual chromosomes with default parameters. The number of potential MQTLs is then estimated on the basis of the model choice criteria from AIC, the modified AIC (AICc and AIC3), Bayesian information criterion (BIC), and Approximate Weight of Evidence (AWE). A model with the minimum values of the selection criteria in at least 3 out of the 5 models was selected as the optimal MQTL model. In the next stage, the 95% CIs and positions of each MQTL were obtained as per the optimal model selected in the first stage. The QTLs were integrated in a way that the peak position of the initial QTLs lay in the MQTL confidence interval. The QTLs whose probability of membership in an MQTL exceeded 60% were assigned to the same MQTL.
Moreover, the QTL-overview statistic proposed by Chardon et al. (2004) was calculated to measure the contribution of a chromosome region to trait variations. The "pnorm" function in R software (R Core Team 2018) was used to obtain the QTL-overview statistic through a step-by-step calculation of the uniform probability that a 0.5-cM-long segment (x and x + 0.5) included a QTL in a test as per Eq. (5).
nbE in which nbE represents the total number of experiments, S i 2 the variance position of the individual QTLs on a chromosome, and nbQTL the number of QTLs. The value of the QTL-overview was obtained in a given chromosome interval based on (i) the number of QTLs lying in the vicinity of the chromosome interval and (ii) the power and accuracy in mapping for these QTLs (high R 2 value or small CI). To observe genomic zones with a significant QTL peak, the mean statistic [U(x)] and a threshold for high values [H(x)] were empirically calculated as 5 times of the mean value . Equations (6)-(7) were, respectively, used to calculate U(x) and H(x).

Identifying the genes located in the MQTL regions
To find the genes located in the identified MQTLs, the related flanking markers were detected on the Oryza sativa genetic map (Temnykh et al. 2001) (http:// archi ve. grame ne. org/ marke rs/ micro sat/) (Supplementary Table S2). The physical positions were then determined after mapping the flanking markers onto the Oryza sativa Japonica group (IRGSP-1.0) reference genome (Kawahara et al. 2013). The genes located in the MQTL regions were ultimately detected in biomart on the ensemble website (https:// plants. ensem bl. org/ bioma rt/ martv iew/) (Supplementary Table S3).

Graphical representation
The distribution and position of the MQTLs on the individual chromosomes were presented as a heat-map based on the rice genome by the ggplot2 R package (Wickham et al. 2019). Moreover, a graphical summary of initial QTLs, MQTLs and QTL-overview statistics was drawn on all the 12 rice chromosomes by SOFIA package (Diaz-Garcia et al. 2017) in R environment.

Gene ontology enrichment analysis
The genes located in the MQTL regions were functionally classified using the web-based AgriGO 2.0 (systemsbiology. cau.edu.cn/agriGOv2/) by the singular enrichment analysis (SEA) tool based on the setting of parameters as follows: (i) p value < 0.05 as the level of statistical significance and (ii) the Fisher's exact test using the adjustment method proposed by Benjamini-Yekutieli for controlling the false discovery rate in multiple tests under dependency. Moreover, the singular enrichment analysis was performed to detect gene ontology terms, i.e. molecular functions, biological processes and cellular components, which were significantly enriched by genes for the individual traits. Gene ontology (GO) enrichment analyses were performed for the genes located in the 64 identified RSA associated MQTL regions, 36 MQTL regions with an interval of less than 1 Mb, the differentially expressed candidate genes (DECGs), and constitutively expressed candidate genes (CECGs), independently.

Collection of datasets and analysis of gene expression
Several DEGs were collected from different microarray (5 published articles) and RNA-seq (15 published articles) experiments at https:// www. ncbi. nlm. nih. gov. The genes with a cutoff of log2-fold of over 1 (two-fold absolute value) and a P-value of at most 0.05 were considered differentially expressed between the genotypes (Supplementary Table S4 and S5). Venn diagram was used to compare the rice root drought-responsive genes detected by RNA-seq or microarray data analysis, and the genes located in MQTL regions. Candidate genes shared between MQTL, microarray and RNA-seq data were analyzed in MapMan . Four overviews, i.e. regulation, metabolism, proteasome and transcription, were used to describe any genes up-regulated in response to drought stress.

Collecting RSA-related genome-wide association studies (GWAS) and comparing them with MQTLs
GWAS studies for root morphological traits in rice were reviewed (Bettembourg et al. 2017;Biscarini et al. 2016;Courtois et al. 2013;Kadam et al. 2017;Li et al. 2017;Mai et al. 2020;Pariasca-Tanaka et al. 2020;Phung et al. 2016;Wang et al. 2018;Xu et al. 2020), and the reported SNP peak positions were collected to find the overlaps between their positions with MQTLs. The genes located in the SNP peak positions (± 25 kb) were extracted from Oryza sativa Japonica group (IRGSP-1.0) reference genome according to their physical positions.

Genetic consensus map construction
A novel integrated consensus genetic map comprising 5447 markers was constructed in BioMercator V4.2 using 1923 SSRs, 3223 RFLPs, 25 ESTs, and 152 other genomic loci based on 5 previously published genetic maps (Supplementary Table S1). The total length of this consensus map was 1528.3 cM, with an average chromosome length of 127.4 cM and a range of 83.0 cM (Chromosome 10) to 181.8 cM (Chromosome 1). The mean number of markers in individual chromosomes was obtained as 454, with the lowest and highest numbers of markers being, respectively, associated with chromosome 10 (n = 256) and chromosome 1 (n = 833). The marker density was obtained as 2.6 to 4.6 markers per cM on chromosomes 12 and 1, respectively, with an average of 3.6 per cM.

Distribution of initial QTLs associated with RSA traits in rice
The reports in the literature and rice database (http:// qtaro. abr. affrc. go. jp) were reviewed to collect QTLs data on the RSA of rice in normal and drought conditions (  Table S6). One hundred and forty-three QTLs were not used in the meta-analysis because either there was no common marker between individual genetic maps and the consensus genetic map, or some of QTLs showed a large confidence interval. Among the 425 initial QTLs, 233 (54.82%) and 192 (49%) were found under normal and water deficit conditions, respectively (Fig. 1a). The present study investigated different populations of QTL mapping studies on root traits, which included F 2 (2 populations), BC (11 populations), RILs (17 populations), DH (5 populations) and NILs (8 populations) ( Table 1). The number of markers used in the previous reports of QTL mapping studies ranged from 215 (Courtois et al. 2003) to 584 (Catolos et al. 2017). The distribution pattern of the initial QTLs on all the twelve rice chromosomes showed that the highest frequency of initial QTLs was related to chromosomes 1 (n = 60), 4 (n = 59) and 7 (n = 58) and the lowest to chromosomes 10 (n = 5), 12 (n = 6) and 5 (n = 16) ( Table 2 and Fig. 1b). The QTLs were unevenly distributed on rice chromosomes, with different combinations of QTLs for different root traits. The average number of initial QTLs for each individual RSA trait was 39, ranging from 13 for the root-to-shoot ratio to 97 for the root length, followed by 65 (15.3%) initial QTLs for root weight traits (i.e. root fresh weight and root dry weight) ( Fig. 1b and Table 3). With an average value of 14.83 cM, the 95% confidence intervals (CI) varied between 1.18 and 58.40 cM, of which approximately 42% of the collected initial QTLs had a CI lower than 10 cM, and 75% had a CI lower than 20 cM (Fig. 2a).
The proportion of phenotypic variance explained (PVE) by the single QTLs ranged from 2.0 to 66.6% with an average of 13.4% (Fig. 2b). The proportion of phenotypic variance explained by each of the 425 initial QTLs was used to rank them in terms of all the RSA traits ( Fig. 2c). A total of 271(59.3%) of the 425 QTLs showed a PVE of more than 10%, whereas 186 QTLs (40.7%) explained less than 10% of the phenotypic variance. For the root length, among the 97 represented QTLs, 87 QTLs explained more than 10% of the phenotypic variance ( Fig. 2c).

Map projection and meta-analysis of QTLs controlling RSA traits
The collected initial QTLs were projected on the integrated consensus map with its high-density SSR and RFLP markers. A total of 64 MQTLs were developed from the 425 initial QTLs for RSA traits with at least 2 MQTLs on each of the twelve rice chromosomes. The number of MQTLs per chromosome ranged from 9 MQTL on chromosome 2, to 2 MQTL on chromosomes 10 and 12 (Figs. 2d, 3 and 4; Table 3). Integrated QTL data showed that all rice chromosomes seemed to be involved in genetic control of root traits. The frequency of clustered initial QTLs per MQTL lay between 2 QTLs in 8 MQTLs (MQTL1-8, MQTL2-1, MQTL6-2, MQTL7-4, MQTL9-2, MQTL10-1, MQTL11-2 and MQTL12-1) and 22 QTLs in MQTL3-1 on chromosome 3 (Table 3). Forty-four (67.8%) MQTLs were obtained by clustering QTLs from at least 3 different experiments, which involved diverse types of the mapping population. The stability of these MQTLs was more likely to be preserved in a variety of environments. The mean of phenotypic variance explained (PVE) by each MQTL ranged from 6.8 to 30.7% with an average of 13.6%, while the 95% CIs reported for the MQTLs ranged from 0.11 cM (0.16 Mb) for the intervals C847-RM473A on chromosome 7 to 14.23 cM (7.98 Mb) for the intervals RM6179-RM4455 on chromosome 10 ( Table 3).
The confidence interval was narrower in the individual MQTL regions than the mean CI of the original QTLs in that position. At 17 MQTL regions, the CI was reduced to less than 2 cM, with a reduction in length by 12.6 times of the mean CI initial QTLs. Moreover, the maximum decrease was observed in the CI of the initial QTLs on chromosome 7, whereas the MQTL7-7 CI was by 136.9 times lower than the CI mean in the clustered QTLs and corresponded to a confidence interval of merely 0.11 cM (0.16 Mb). A total of 14 MQTLs had the physical length of around 500 kb and 10 of which also had less than 2 cM genetic distance (Table 3). These ten MQTL regions also had a mean phenotypic variance of 14.8%. The flanking markers of these ten MQTL regions were appropriate for molecular breeding and marker-assisted selection in future genetic improvement programs of root morphological traits in rice. The physical intervals of 36 (56.2%) MQTLs were below 1 Mb. The physical length of MQTL ranged from 0.08 Mb (MQTL3-5) to 7.98 Mb (MQTL10-1). The number of RSA traits in each MQTL ranged from one (root length) in 1 MQTL (MQTL1-8) to eight in 4 MQTLs (MQTL1-6, MQTL4-3, MQTL4-4 and MQTL7-7). The individual QTLs for the root length and root thickness were present in 49 MQTLs out of the 64 detected MQTL regions, the most for any root morphological traits, while the individual QTLs for the root surface area was only present in 11 MQTLs, the least among the RSA traits (Table 2).

Estimation of QTL-overview index for RSA QTLs in rice
After projecting the 425 QTLs onto the consensus genetic map, the density of the presence of QTLs described as "QTL-overview index" was calculated for the considered interval of 0.5 cM on each chromosome to identify genomic regions significantly associated with RSA traits ( Fig. 3 and Supplementary Fig. S1). Sixty-four out of 73 overview index peaks obtained were higher than 0.0122 as the mean value of the statistic over the genome and showed 'real QTLs' that affect RSA traits in rice. According to Supplementary Fig. S1, fifteen out of the 64 peaks also exceeded 0.061 as the high-value threshold. The number of significant peaks ranged from one in chromosomes 1, 2 and 5 to three in chromosomes 4 and 9. In MQTLs on the twelve chromosomes in rice. DRR: deep root length ratio, RDR: ratio of deep rooting, RDW: root dry weight, RFW: root fresh weight, RGR: root growth rate, RL: root length, RN: root number, RSA: root surface area, RSR: root to shoot ratio, RT: root thickness, RV: root volume 1 3 the latter case, all the 5 peaks exceeded the high-value threshold. The results revealed that the estimated CIs identified through the QTL-overview analysis were similar to the meta-analysis method. The QTL-overview index is a statistical approach that has been widely used for traits such as flowering time in maize , the response of leaf growth to water deficit in maize (Welcker et al. 2007), yield and related traits in maize (Martinez et al. 2016), the fatty acid content in soybean (Qin et al. 2018) and root-related traits in bread wheat (Soriano and Alvaro 2019).

Detecting differentially expressed genes (DEGs) in the rice root under drought conditions
The drought-responsive genes in the rice root were collected from RNA-seq and microarray datasets (Supplementary Table S5). Based on the RNA-seq data analysis, 10,239 up-regulated and 7214 down-regulated rice root DEGs were observed at drought compared to normal conditions. Microarray meta-analysis suggested 22,307 DEGs, among which, 16,290 and 6017 were up-and down-regulated, respectively (Supplementary Table S7 and Fig. S2). initial QTLs on the twelve chromosomes of rice (blue lines) and position of MQTLs with the 95% confidence intervals (red). d QTL-overview index (probability density) for the root traits QTLs on consensus genetic map. e Proportion of phenotypic variance explained (R 2 ) for each initial QTLs. f Density of initial QTL for each MQTL Finally, MQTL regions associated with RSA were explored to find the drought-responsive genes in the rice root. A total of 5448 and 1338 common genes were revealed using the Venn diagram between the DEGs derived from RNA-seq and microarray data and the genes located in MQTL regions (Supplementary Table S8 and Fig. S2) for all the 64 MQTLs ( Supplementary Fig. S2b) and the 36 MQTLs with CI of less than 1 Mb ( Supplementary Fig. S2a), respectively. These genes that are both differentially expressed at drought conditions in rice roots and located in MQTL regions can be considered as differentially expressed candidate genes (DECGs).

Gene ontology analysis of the genes located in the RSA-associated MQTL regions
Enrichment analysis of the genes located in the 64 MQTLs, 36 MQTLs, DECGs and CECGs, independently revealed many shared and a few diverged enriched GO terms (Supplementary Table S9 and Fig. S3). It confirms the significance of the results and provides evidence that genes with associated functions might be clustered and contributed to the QTL traits. The common terms included gene expression, protein modification process, RNA metabolic process, phosphotransferase activity, protein serine/threonine kinase activity, c Distribution of initial QTLs on twelve chromosomes of rice (black lines) and position of MQTLs with confidence intervals (red areas). d Coefficient of reduction in physical CI from mean original QTLs to MQTL. (e) Number of genes involved in each MQTL interval calcium ion binding transporter activity, nucleus, membrane part and ribosome. However, as expected, some terms were enriched in the drought-responsive candidate genes, but not in the constitutively expressed candidate genes, such as ARF protein signal transduction, response to stress, ion transport, response to oxidative stress, and plant-type cell wall organization. Under the gene expression term, many RSA-related genes were observed such as members from WRKY, NAC, ARF, WOX, AUX/IAA, RR, RAA, and HSFA2 families, which were also detected by MapMan as drought-responsive TFs. Likewise, other significant GO terms such as the protein modification processes were significant through both Map-Man and GO analysis and contained RSA-related genes (e.g. ACA8, SAPK 6,10,7 and CCR3).

Analyzing the drought-responsive candidate genes in MapMan
The fold-change data and Locus IDs for 5448 and 1338 drought-responsive genes of the 64 and 36 MQTLs (with an interval of less than 1 Mb), respectively, were uploaded to MapMan toolkit (Supplementary Table S10) and the various overviews were obtained ( Supplementary Fig. S4a to S4f). It is worth noting that almost similar results were obtained from MapMan analysis of the DEGs located in the 64 and 36 MQTL regions. Investigating the regulation overview of the 5448 genes showed the up-regulation of 417 transcription factors (TFs), 276 genes related to protein degradation, and 192 genes associated with protein modification in the rice root under drought stress ( Supplementary Fig. S4b). The highest frequency of TFs and their importance in the regulation of the root response and structure under drought stress, warrant further investigations ( Supplementary Fig. S4c).
The most significantly enriched genes in the signaling pathway included those for G-protein and calcium regulation. Thioredoxin as the reduction-oxidation response element dominated the cellular processes of drought-responsive genes in rice. A G protein-mediating calcium-signaling cascade, therefore, functions upstream of drought-responsive TFs in the hierarchy of signaling pathways to regulate the reduction-oxidation reaction through the thioredoxin activity and trigger the drought tolerance. This evidence suggests the cooperative role of thioredoxin, G-proteins and calcium regulation in responding to drought in the rice root. It is worth noting that MapMan analysis results of the DEGs located in the 64 and 36 MQTL regions were somehow similar ( Supplementary Fig. S4a to S4f).

Validation of MQTLs with the related GWAS studies
Significant overlaps existed between the MQTLs detected using meta-analysis and the SNPs detected using the GWAS method that are linked to RSA traits in rice genome.
Interestingly, 52 MQTLs out of the 64 identified MQTLs were co-located with 171 SNP peak positions reported in rice GWAS for root morphological traits (Figs. 5 and 6). A total of 752 rice genes were located in the SNP peak positions (± 25 kb) overlapping with MQTLs, among them 49 SNPs were occurred exactly within a gene (Supplementary  Table S11).

MQTLs for RSA traits and the perspective of their implications in MQTL-assisted breeding
The RSA plays a key role in absorbing nutrients and water by a plant, its anchoring to the soil and its yield (Subira et al. 2016). Linkage analysis and linkage disequilibrium mapping indicate that root morphology in rice is a complex trait that involves multiple loci with small effects. Metaanalysis is a powerful statistical technique that can be performed on QTL data collected from different genetic backgrounds and environments to determine consistent QTLs, enhance the accuracy of their chromosomal locations, and integrate the QTL data collected from numerous independent experiments (Goffinet and Gerber 2000). In the present QTL meta-analysis, we collected the data of 425 QTLs from 34 articles published between 2001 and 2020 (Supplementary Table S6) to identify genomic regions linked to RSA and the drought stress tolerance in rice. We could detect 64 MQTLs by performing a meta-analysis based on the modified Akaike information criterion (AIC) ( Table 3). The supporting intervals of the identified MQTLs with an average of 3.79 cM (1.57 Mb) reduced 3.88 times as compared to the mean of corresponding original QTLs. The only similar previous work (Courtois et al. 2009) just covered 306 QTLs reported until 2007, detecting 119 MQTLs with 4.14 Mb average CI for the initial QTLs and 1.98 Mb for MQTLs. In our study, 376 out of the 425 initial QTLs were SSR flanking markers, of which 268 QTLs belonged to the reports after 2007. In other QTL meta-analyses on diverse quantitative traits in various crops, including bread wheat, barley, maize and soybean, has been reported a reduction of 10%-21% in the total number of initial QTLs compared to the number of MQTLs, and the average reduction in the CI of the MQTLs varied from two to four times of the initial QTLs (Ballini et al. 2008;Courtois et al. 2009;Darzi-Ramandi et al. 2017;Hao et al. 2010;Khahani et al. 2020;Lanaud et al. 2009;Rong et al. 2007;Soriano and Alvaro 2019). In the current research, the genetic CI of 71.8% and physical CI of 75% MQTL's were narrower than 5 cM and 2 Mb, respectively. The present study found several MQTLs in recombination hot spot regions with a narrow CI such as MQTL7-7 (0.16 Mb, 0.11 cM), MQTL3-5 (0.08 Mb, 1 3 0.22 cM), MQTL6-6 (0.27 Mb, 0.47 cM), MQTL9-6 (0.27 Mb, 0.48 cM), MQTL5-3 (0.21 Mb, 0.71 cM) and MQTL7-2 (0.51 Mb, 0.82 cM) to be useful for the finemapping studies of quantitative trait loci (Table 3). The present findings suggested that analyzing MQTLs can help discover consensus genomic regions and accurately locate linked QTLs for RSA traits in rice.
Although there is still a high diversity for the quantitative inheritance nature of RSA traits, rice breeders can take advantage of the current genetic diversity to select the most appropriate MQTLs for improving the root characteristics and drought tolerance in their elite breeding materials. According to Löffler et al. (2009), the MQTLs selected for breeding purposes should have (i) a narrow CI, (ii) comprise a large number of initial QTLs, and (iii) a high proportion of its mean phenotypic variance is explained by initial QTLs. Using some of the identified MQTLs for MQTL-assisted breeding to improve the RSA traits to enhance drought Bettembourg et al. 2017 Biscarini  Table S11) tolerance appears promising. Interestingly, eleven of the identified MQTLs contained at least 10 initial QTLs with a 95% CI of below 2 cM (Table 3). Among them, the MQTL1-7, MQTL3-5, MQTL4-3, MQTL4-5, MQTL6-6, MQTL7-7, and MQTL9-6 span less than 500 kb and the original QTLs constitute six to nine independent studies. These MQTLs could be used as multi-effect hotspot regions affecting RSA traits, which is worth applying in breeding programs in the near future.

Potential candidate genes involved in RSA and drought tolerance
Multiple genes associated with root-related traits in rice were found to be located in MQTL intervals (Supplementary Table S3) which were used to reconstruct the related molecular network, and will be discussed in the following ( Fig. 7 and Supplementary Table S4).

Plant hormone-mediated root growth and development
We observed several genes involved in auxin biosynthesis or signaling located in the MQTL regions ( Fig. 7 and Supplementary Table S4). Auxins regulate many growth and development dimensions of root through elongating primary roots and root hairs and raising the number of lateral root primordia (Overvoorde et al. 2010). Some members of the rice YUCCA gene family involved in the biosynthesis of indole-3-acetic acid (IAA), including OsYUCCA1 located in MQTL1-5, OsCOW1, OsYUCCA8 and OsNAL7 located in MQTL3-1, play a key role in the root growth (Fujino et al. 2008;Woo et al. 2007). It is reported that overexpression of a YUCCA gene, which encodes the rate limiting enzyme of auxin biosynthesis, significantly increases the proliferation of crown roots depending on the presence of WOX TFs . The WOX6 transcription factor lays in MQTL3-3. The WOX transcription is activated by auxin, which then initiates the crown root development in rice by establishing the YUC-Auxin-WOX module . WOX was also found to interact with ERF and bind to the OsRR2 (MQTL2-3) promoter as a response regulator of cytokinin signaling to regulate crown root morphological traits, root number, root to shoot ratio and root volume . As a RING-H2 membrane-anchor E3 ubiquitin ligase, EL5 (MQTL2-5) preserves cell viability after initiating root primordia through cytokinin-mediated signaling (Koiwai et al. 2007).
A rice zinc finger protein (OsZFP) in MQTL1-2 interacts with a cyclophilin (LRT2/OsCYP2) and influences root development through the IAA pathway . LRT2 (LATERAL ROOTLESS2), the cyclophilin protein (OsCYP2) located in MQTL2.1 regulates auxin signaling and lateral root initiation in rice ). LRT2/OsCYP2 was found to contribute to degrading AUX/ IAA proteins (Kang et al. 2013). Several OsAUX/IAA15 and OsAUX/IAA22 members are found in MQTL5-1 and MQTL6-3 regions, which act as the negative regulators of ARFs (Mockaitis and Estelle 2008). Based on our GO analysis results, regulation of ARF GTPase activity and ARF protein signal transduction were enriched in the biological processes. The activity of ARF GTPase activator was also enriched among the molecular functions. Two major protein families that regulate auxin signaling include ARF proteins as the DNA-binding transcriptional regulators of auxin responses and, Aux/IAA proteins as the negative regulators of ARF (Mockaitis and Estelle 2008). The rice lateral root formation requires Auxin/IAA protein and ARF-mediated signaling (Yamauchi et al. 2019). OsARF1 was found to lie in MQTL5-2, OsARF7 in MQTL2-5 and OsARF16 in MQTL6-2. It is reported that OsARF7 and OsARF16 regulate RSA through AUX/LAX in auxin signaling (Lee et al. 2019;Mockaitis and Estelle 2008;Yamauchi et al. 2019). As functional auxin influx carriers, the AUX/LAX family of proteins mediate auxin-related developmental programs. Responses to abiotic stress, root gravitropism and root hair development are controlled by AUX1. Lateral root development is affected by LAX3 and AUX1 (Swarup and Bhosale 2019;Swarup and Péret 2012). OsAUX3 is expressed in root hairs, primary and lateral roots. OsAUX3 mutations were found to decrease the length of primary roots and lateral root density and increase the length of root hairs (Wang et al. 2019a). OsAUX3/LAX2 was located in MQTL3-2. Adventitious root formation in Arabidopsis to be affected by LBD18 and LBD16 functioning downstream of ARF19 and ARF7 (Lee et al. 2019). We found LBD38 gene in the MQTL3-4 region that is a class II type LBD protein, acting as a transcriptional activator (Pan et al. 2017).
Oryza sativa root architecture associated 1 (OsRAA1) located in MQTL1-2 play a role in the auxin-mediated development of root in rice (Ge et al. 2004). The transition from metaphase to anaphase during cell division is inhibited by OsRAA1, and its degradation by the ubiquitin-proteasome system is essential for initiating anaphase in mitosis during root development (Xu et al. 2010). OsMADS57, a MADSbox transcription factor, was found in MQTL2-8. The inhibition of seminal root elongation in OsMADS57 mutants can be explained by elevated auxin and polar auxin transport to the root tips of the mutant plants (Huang et al. 2019). MAIF1, a rice F-box domain gene located in MQTL2-5, was found to contribute to several signaling pathways in the regulation of root growth. Auxin, abscisic acid, cytokinin and abiotic stress can induce the expression of MAIF1. The root growth in rice can be stimulated through the over-expression of MAIF1 (Yan et al. 2011).
OsEIL1, located in the MQTL3-3 interval, activates the OsYUCCA8 transcription and, therefore, auxin biosynthesis, regulating ethylene-inhibited root elongation in rice seedlings (Qin et al. 2017). ACC (1-aminocyclopropane-1-carboxylic acid) oxidase gene of OsACO1 (MQTL9-4) was  Fig. 7 The molecular network associated with rice root development using potential candidate genes located in MQTL regions. An arrow shows a positive regulatory action, and lines with a flat head in their end show negative regulatory actions. Text color code: biological pro-cesses: purple, hormones: blue, genes/proteins located in MQTLs: black, genes/proteins located in both MQTL regions and SNP peak positions reported in rice GWAS for root morphological traits: dark blue, Ca +2 : red, root traits: pink, root parts: green enhanced in transgenic plants overexpressing OsEIL1 (Mao et al. 2006). Ethylene signaling is regulated in rice through activating OsEIL1 by RAVL1 (MQTL4-4), as an upstream component of brassinosteroid biosynthesis and signaling (Zhu et al. 2018). Mediating the degradation of OsEIL1 through the ubiquitination pathway by OsEBF1 (MQTL6-4) as an E3 ligase suggested the negative regulation of the ET signaling pathway in response to the infestation of BPH (Ma et al. 2020). Abscisic acid (ABA) plays a key role in raising the length of root hairs in plants. It can accumulate auxin along the root hair by promoting the biosynthesis and polar transport of auxin in the root tip (Qin et al. 2017). Expression of some ABA biosynthesis genes could regulate early signaling genes, e.g. OsPP2C68 (MQTL9-2), SAPK6 (MQTL2-5) and late responsive gene OsLEA3 (MQTL6-3) during droughts (Xiong et al. 2014). SAPK6 is effective in root elongation during drought conditions (Xiong et al. 2014), and longer root hairs were reported in transgenic rice overexpressing SAPK10 (MQTL3-4) (Qin et al. 2017). Jasmonate (JA) signaling components such as OsJAZ1 (MQTL4-6) was also found, of which OsJAZ1 is reported to play a negative role in regulating drought tolerance through ABA and JA pathways (Fu et al. 2017). The degradation of jasmonate ZIM-domain (JAZ) proteins releases TFs known to interact with JAZ proteins, which can generate responses to growth and stress (Tian et al. 2019). OsJAZ1 was found to be associated with the root weight, and potentially regulate rice root development at different developmental steps (Pan et al. 2017). OsCKI1, a rice casein kinase I detected in MQTL2-6, was found to contribute to different hormone-signaling pathways and regulate the development of rice lateral roots ).

Root growth angle
Deeper rooting 1 (DRO1) (MQTL9-4), regulated through ARF by auxin, appears to contribute to elongating root tip cells causing the asymmetric growth of root and its downward bending in response to gravity. Induction of the DRO1 expression raises the angle of root growth, increasing the downward growth of the root (Uga et al. 2013). Controlling RSA by DRO1 raises the rice yield during drought stress (Uga et al. 2013). Moreover, DRO1 was found to be expressed in the root tip in the vicinity of the root apical meristem as well as in basal shoots in the crown root primordia (Uga et al. 2013). Interestingly, we further found some other candidate genes which might regulate branching and angle in lateral roots; for example, OsLAZY1 (MQTL6-1), OsSOL1 (suppressor of lazy1) (MQTL6-1) and OsTAC3 (MQTL3-5). DRO1-related genes, including TAC1 and LAZY1, belong to the IGT gene family, contain a C-terminal EAR-like motif IVLEI and are observed in different plant phyla (Ashraf et al. 2019). The gravitropism of the root and shoots is affected by LAZY1 and its orthologues in rice, maize, Medicago and Arabidopsis (Ashraf et al. 2019). HSFA2D located in MQTL3.1 functions as a positive regulator of the LAZY1-dependent asymmetric distribution of auxin in the upstream. Auxin causes the asymmetrical expression of WOX6 (MQTL3.3) to connect gravitropism to controlling the tiller angle in rice . PROG1 (PROSTRATE GROWTH1) (MQTL7.1), and LPA1 (LOOSE PLANT ARCHITECTURE1) (MQTL3.2) were reported as the main genetic regulators of the tiller angle in rice (Li et al. 2007).

Lateral root development
The enhancement of lateral root formation causes a potentially useful adaptation to drought in lowland rice under drought stress (Hazman and Brown 2018). A panel of lateral root development-related potential candidate genes was located in MQTL regions ( Fig. 7 and Supplementary  Table S4). O. sativa drought stress response-1 (OsDSR-1), a calmodulin-like gene, was found in MQTL10-1. The binding of its protein to Ca 2+ leads to conformational changes, and increases drought tolerance by scavenging the reactive oxygen species in rice . The expression of OsDSR-1 in Arabidopsis was found to enhance developing lateral root in high K + concentrations and improve sensitivity to ABA (Yin et al. 2011). OsWRKY31 found in MQTL6-3 can act in the auxin signal transduction pathway, and its ectopic expression reduced the length and formation of lateral roots . OsWRKY76 lying in MQTL9-3 was found to regulate cellular responses to abiotic and biotic types of stress (Yokotani et al. 2013) partly by regulating the expression of RSOsPR10 in the lateral roots of rice, improving the mass and growth of root and increasing tolerance to high salt environments and soil desiccation (Yamamoto et al. 2018).

Root diameter
As a transcriptional activator and stress-responsive protein lying in MQTL11-2, OsNAC5 improves stress tolerance through the up-regulation of stress-inducible genes of rice, e.g. DjC10 (MQTL1-4) and OsLEA3 (MQTL6-3). The overexpression of OsNAC5 increases the root diameter in rice and improves grain yield and drought tolerance (Jeong et al. 2013). The root-specific overexpression of OsNAC10 (MQTL7-6) significantly improves the rice yield, especially during drought stress (Redillas et al. 2012). OsNAC10 was found to induce 34 root-specific target genes such as those encoding cytochrome P450, mitogen-activated protein kinase kinases (MAPKK), LEA and TFs such as WRKYs and NACs (Jeong et al. 2010). As a key factor in preserving 1 3 calcium homoeostasis, calcium-transporting ATPase serves as an early response to drought, salinity stress and low temperature in plant cells (Knight 1999). A member of this family, ACA8 (Ca 2+ P-Type ATPase 8) (Knight 1999), was located in MQTL10-2. ACA8 was up-regulated in OsNAC5 overexpressed plants (Jeong et al. 2013).

Multi-level data integration to discover novel candidate genes
A large number of studies have shown that QTL mapping and genome-wide association study (GWAS) are two complementary methods for mapping causal genes and dissecting the genetic basis of the traits of interest (Mahuku et al. 2016;Tao et al. 2013). Meta-analysis is an effective approach to integrate QTLs from several independent experiments, and GWAS is a robust method for detecting significant effects in diverse germplasms, each has its own merits and demerits, which can complement each other (Sonah et al. 2015). A major challenge of QTL analysis is the reduction of the CI to discover the most relevant genomic regions, and the integration of GWAS results with meta-analysis of QTLs could help to do so very effectively. In our study, significant overlaps existed between the MQTLs detected using metaanalysis and the SNPs detected using the GWAS method linked to RSA traits in rice genome. Out of the 64 identified MQTLs, 52 MQTLs were overlapped with 171 SNP peak positions reported for root morphological traits (Figs. 5 and 6). Among them, all MQTLs on chromosomes 1, 2, 8, 10 and 12 were in co-linear region with GWAS signals. This integrative approach leads to the detection of 49 rice genes located in both the MQTL regions and SNP peak positions reported in rice GWAS for root morphological traits including genes encoding EXPA6, MADS59 and RR3. In addition, looking at a 25 kb region on either side of the SNP peak positions, 755 genes were found comprising many novel genes and some genes with known RSA-related functions such as DRO1, WRKY62, CKI1, ARF8, SAUR30 and HOX11 (Supplementary Table S11).
To find the chromosome regions with the most contribution to a desired phenotype, the QTL-overview index is calculated based on the number of original QTLs, accuracy in mapping for these QTLs and the high phenotypic variance ). In the current research, the genes located in the significant peaks of both SNPs and the QTL-overview for root architecture traits were considered as novel candidate genes (Table 4). These 10 novel genes are interesting candidates for functional analysis, while the transcript level of rNBS41, GELP58, EME1 and bZIP64 was also drought-responsive in the rice roots.

Conclusion
Root system architecture plays critical roles in plant growth and productivity as well as plant tolerance to abiotic stresses such as drought. Therefore, improving RSA through molecular breeding or genetic engineering is a promising strategy to develop the yield and stress tolerance in rice. In the present research, a genome-wide meta-analysis of RSA in rice led to the identification of 64 MQTLs. Interestingly, seven of these MQTLs have a physical length of around 500 kb and genetic distance of less than two cM, contained at least ten initial QTLs while the original QTLs were derived from six to nine independent studies, and cover a mean phenotypic variance of 14.8%. Two of them (MQTL1-7 and MQTL6-6) were further confirmed with SNP peak positions reported Table 4 The novel candidate genes located within both the SNP peak positions and the QTL-overview peaks for root architecture traits Trait Genes (RAP ID) Gene name Gene annotation Gene start (bp) Gene end (bp) Chr SNP position ( 1 3 in rice GWAS for root morphological traits. Conclusively, these MQTL regions can be appropriate for MQTL-assisted breeding in future genetic improvement programs of RSArelated traits in rice. On the other hand, the DEGs in the rice root under drought conditions were detected through the analysis of RNA-seq and microarray datasets, and the MQTL regions associated with RSA were explored to find the drought-responsive genes in the rice root. Gene ontology enrichment analysis of the genes located in these MQTLs provides evidence that the genes with RSA associated functions are clustered and contributed to the QTL traits. Multiple genes with functions associated with RSA and drought tolerance were found in the MQTLs comprising some genes from YUCCA , WOX, RR, LRT2/CYP2, AUX/IAA, ARF, RAA, MAIF1, EIL, PP2C, SAPK, DRO, DSR, WRKY, RSOsPR, NAC, C4H, PAL, PRX and CAD, EXPA, GLU and ARD families. The genes located within both the SNP positions and in the QTL-overview peaks for root architecture traits are suggested as interesting candidate genes for further functional analysis. Conclusively, we integrated meta-analysis of QTLs associated with RSA traits in normal and drought stress conditions, GWAS studies results for root morphological traits, and the transcriptome analysis data of the root under drought conditions leading to find 189 candidate genes (Supplementary Fig. S5), which might play an important role in rice RSA under drought stress. After functional analysis, the promising candidate genes would be applicable for root genetic engineering aimed at improving yield potential, stability and performance in drought conditions.