Distribution of initial QTLs on nine chromosomes in foxtail millet
A total of 448 QTLs were collected from ten published studies on yield in foxtail millet. The studies included ten diverse experimental crosses with 11 parental lines and 2,092 progeny lines and various population types including three F2 and seven recombinant inbred lines (Table 1), while the population size ranged from 124 (Wang et al. 2017) to 543 individual genotypes (Wang et al. 2019). The original QTLs were dispersed irregularly on all the nine chromosomes of foxtail millet (χ2 = 439.3), the highest number of QTLs were found on chromosome 9 with 90 QTLs, followed by chromosomes 5 (n = 80), 7 (n = 58) and 1(n = 57) and chromosome 8 had the fewest QTLs with 19 QTLs.
Among these 448 initial QTLs, 351 (78.3%) and 97 (21.7%) were found under normal and drought stress conditions, respectively (Fig. 1a). Under normal conditions, the number of QTLs varied from 16 QTLs on chromosome 8 to 66 QTLs on chromosome 5. The distribution of QTLs under drought conditions was also diverse on various chromosomes; 25 QTLs were located on chromosome 9, and the lowest number of QTLs (4 QTLs) were located on chromosome 1. The number of QTLs per trait ranged from 7 to 118, where tiller number (TN), plant height (PH), and biomass (BM) exhibited the highest number of initial QTLs with 118, 92, and 62 QTLs, respectively (Fig. 2). Moreover, the distribution of QTLs for each trait varied across environments; TN and PH QTLs were mainly evaluated under normal conditions, with 114 and 80 QTLs, respectively. The number of initial QTLs reported under drought stress conditions ranged from 1 (panicle length) to 46 (water status) for each trait.
The 95% confidence interval (CI) for each QTL differed from 1.45–49.05 cM with a mean of 18.37 cM. Approximately 89 (19.8%) of the collected QTLs had a CI less than 10 cM, whereas CI of 196 QTLs (43.7%) were between 10 to 20 cM (Fig. 1d). The average PVE for the initial QTLs was 7.5%, with a maximum of 61.4% (Fig. 1c). Among the 448 initial QTLs, 26 QTLs (5.8%) exhibited a PVE greater than 15% (Fig. 1b), and the logarithm of odds ratio (LOD score) ranged from 2.03 to 49.75, with a mean value of 5.8 for combined normal and drought stress in foxtail millet.
In normal conditions, the average value of 95% confidence intervals (CI) was 20.65 cM, with approximately 11.68% of the collected original QTLs demonstrating a CI lower than 10 cM, and 53.2% a CI lower than 20 cM. The PVE by the single QTLs fluctuated from 1.99 to 61.4% with a mean of 6.3%. Among the 351 initial QTLs detected under normal conditions, 297 QTLs (84.6%) indicated a PVE less than 10%, while 54 QTLs (15.3%) explained more than 10% of the phenotypic variance, and the logarithm of odds ratio (LOD score) ranged from 2.03–49.75 with a mean value of 5.72 (Table S1). Among the 97 initial QTLs found under drought stress, 48 QTLs (49.4%) exhibited a CI of less than 10 cM, and 85 QTLs (87.6%) exhibited a CI of less than 20 cM. PVE value explained by these initial QTLs ranged between 4.96 and 27.83 (average: 10.86%). The range of LOD scores for drought tolerance QTLs was 2.03 to 19.52, with a mean of 5.18.
Drought-tolerance Associated Meta-analysis Of Qtls
The collected 448 individual QTLs were mapped to the consensus genetic map, and 41 meta-QTL (MQTL) regions were identified on the nine chromosomes of foxtail millet (Table S3). The number of MQTLs per chromosome extended from seven MQTLs on chromosome 2 to two MQTLs on chromosome 8 (Table S3).
On chromosome 9, MQTL9-5 was obtained with a relatively narrow CI (5.07), with 34 initial QTLs reported from three experiments containing the highest number of QTLs, followed by MQTL7-5, MQTL1-4, and MQTL9-3 with 26, 21, and 20 initial QTLs derived from eight, two, and three experiments, respectively, and covering a relatively narrow CI of 3.18, 0.31, 7.95, and 1.13 cM, respectively (Table S3). Inclusively, MQTL5-6 with 19 initial QTLs had the highest mean percentage of phenotypic variation (PVE = 15.79) containing five studies for five traits, including PH, WUS, TN, YLD and RL, followed by MQTL9-3, MQTL5-5, and MQTL9-1 with 14.23%, 12.57%, and 11.14% PVE value, which can be regarded as the most efficient QTLs for the associated traits (Table S3). Notably, these MQTLs support the important traits, for instance, MQTL9-5 containing six studies for nine traits including BM, PH, NN, PL, SD, WUS, TN, YLD and RL, and MQTL7-5 consisting of eight studies for seven traits including BM, PH, NN, PL, SD, TN, WUS (Table S3).
Under normal conditions, 38 MQTLs were found on all the chromosomes (Table S3). MQTL5-1 with 23 original QTLs, had the highest number of QTLs, followed by MQTL9-2 and MQTL1-4, with 22 initial QTLs, and CIs of 2.83, 1.98, and 6.5 cM, respectively (Supplementary Table S9). MQTL5-5, with 14 initial QTLs had the highest PVE (17.30%), containing two studies for three traits including PH, BM, and YLD, followed by MQTL5-4 with 11.85% PVE, which can be considered the most effective QTL for the involved traits. Some MQTLs for different traits were identified under normal conditions; for instance, MQTL5-1 consisted of four studies for six traits, including BM, PH, TN, NN, PL and SD, and MQTL9-2 consisted of three studies for five traits, including YLD, TN, PH, BM, and NN.
Under drought stress, 97 individual QTLs were mapped to the consensus genetic map, and 9 MQTLs were found (Table S3). MQTL7-1, with ten initial QTLs, had the most QTLs, followed by MQTL9-5 and MQTL9-3, with nine and eight initial QTLs, and a relatively narrow CI of 1.37, 3.13, and 0.69 cM, respectively. With eight initial QTLs, MQTL9-3 had the highest PVE (22.6%). The physical intervals varied from 0.40 Mb (MQTL7-5) to 16.46 Mb (MQTL1-3), averaging 3.28 Mb (Table S3). MQTL9-5 contained three distinct traits, namely WUS, RL, and BM, followed by MQTL7-1, which contained four studies for four traits including BM, PL, PH, and WUE (Table S3).
Estimation Of Qtl-overview Index
The QTL-overview index was considered for cM units over all chromosomes (Fig. S2). A total of 33 overview peaks were achieved, of which 20 exceeded the average value of the statistic for each chromosome and represented “real QTL” regions. Additionally, seven of the 33 peaks exceeded the high-value threshold (measured as 0.0743 for a total of nine chromosomes) and represented “QTL hotspot” regions. These seven overview peaks corresponded to seven of the 41 MQTLs (17.03%) and contained 99 initial QTLs (22.09%). Two of the seven MQTLs on chromosome 2 (28.5%) had a high overlap with overview peaks. Moreover, after mapping the 50 QTLs onto the chromosome 7 consensus genetic map, one of the five MQTLs (20%) was detected in hotspot overview peaks. A total of six overview peaks were obtained for 68 QTLs on chromosome 5; two of six (27 QTLs) MQTLs were located in the overview peak interval. Moreover, four overview peaks were identified on the chromosome 9, of which two (33 QTLs) were identified based on meta-QTL results (Fig. S2).
Distinguishing Differentially Expressed Genes (Degs) In Foxtail Millet Under Drought Conditions
The drought-responsive genes were collected from RNA-seq data of six independent studies and one microarray experiment (Table S8). According to the RNA-seq data analysis, 9,356 DEGs (6,016 up-regulated and 3,340 down-regulated) were found under drought stress compared to normal conditions. In addition, the microarray assay suggested 26 down-regulated and 74 up-regulated DEGs (Table S5).
A total of 1,631 and 314 common genes were discovered using the Venn diagram among the DEGs derived from RNA-seq and microarray data and the genes located in the MQTL regions for all the 41 MQTLs and the 10 MQTLs with CI lower than 1 Mb, respectively (Fig. S3; Table S9).
Functional Annotation Of Degs Located In The Mqtl Regions
Among the 8564 genes located in 41 MQTLs, 1,631 were differentially expressed candidate genes and 6,933 constitutively expressed candidate genes (Table S10, Fig. S4). The Singular Enrichment Analysis (SEA) were used to inspect the gene ontology enrichment using the DEGs list located in 41 identified MQTL regions (1,631 DEGs out of 8564 genes) and 10 MQTL regions with an interval of less than 1 Mb (314 DEGs out of 1,605 genes).
The results of the GO analysis showed that 1,070 of the 1,631 genes were functionally annotated using the AgriGO v.2.0 database (Table S11). The number of annotated genes ranged from 0 to 73 (MQTL5-4), and the percentage of significant GO terms varied from 0 to 65.12. Only 19 out of 40 MQTLs contained a total of 231 significant GO terms. In the biological process category (BP), GO analysis of the DEGs found in the MQTL regions shown that signaling, metabolic process, biological regulation, cellular process, and response to stimulus were significantly enriched (Fig. S4). Several molecular function categories were significantly enriched, including binding, transcription regulator, structural molecule, electron carrier, catalytic, and transporter activity. The significantly enriched terms of cellular component (CC) category, were comprised membrane-enclosed lumen, organelle, symplast, macromolecular complex and extracellular region (Fig. S4).
Metabolic Pathways Of Drought-regulated Genes
To supplement the GO analysis, we used MapMan analysis to search for the putative functions of the identified drought-response genes in the MQTL regions. The MapMan file contained the fold-change data and locus identifiers for 1,631 drought-responsive genes located in the 41 MQTLs (Table S12). The overview of regulatory pathways of the 1,631 DEGs indicated on the up-regulation of 86 transcription factors (TFs), 27 genes associated with protein degradation, and 48 genes related to protein modification in foxtail millet under drought stress (Fig. S5a). Genes involved in hormone signaling pathways including IAA, ABA, ethylene, and jasmonic acid, were also enriched (Fig. S5a). Additionally, the regulation overview revealed the identification of 17 genes for calcium regulation, 44 genes for receptor kinase, six genes for G-proteins, one gene for phosphoinositide, and five genes for light signaling pathways (Fig. S5a).
The overview of stress response pathways showed that genes involved in abiotic stress response, signaling, redox state, peroxidases, and glutathione S-transferases were enriched (Fig. S5b). Genes involved in cell wall modification (xyloglucosyl transferase, xyloglucan endotransglycosylase) and cell wall degradation (polygalacturonase) represented another set of important components that were enriched (Fig. S5b). The results of mapping the DEGs to secondary metabolite pathways overview revealed that lignin, terpenoids, flavonoids, and phenylpropanoids pathways were among the enriched pathways (Fig. S5c). The cellular overview pathway indicated that genes encoding abiotic stress-related miscellaneous enzyme families (misc) and development were up-regulated in foxtail millet under drought stress (Fig. S5d). Although the genes involved in transcription pathways were mapped in 41 MQTLs, our results uncovered that most of the TF families belonged to the eight myeloblastosis (MYB), nine Apetala2/Ethylene Responsive Element Binding Proteins (AP2/EREBPs), nine Basic Helix-Loop-Helix (bHLH) genes, seven homeobox genes (HB), six Chromatin remodeling factors, two SET-domains, 11 WRKY domains, three MADS-box genes, 10 Histone genes, 7 C2H2 zinc finger family, five putative DNA-binding, 4 G2-like transcription factor family (GARP) and four basic leucine zipper (bZIP) (Fig. S5e).
Ortho-mqtls In The Genomes Of Barley, Maize, Rice, And Foxtail Millet
In order to investigate ortho-MQTLs, 41 MQTLs derived from foxtail millet were compared to the genomic regions related to drought tolerance and yield associated traits reported in barley, maize, and rice. Consequently, four, two, and four ortho-MQTLs were identified for foxtail millet and barley, maize, and rice, respectively (Table 2, Fig. 4). As a result, foxtail millet MQTL2-5, MQTL3-1, and MQTL5-7 were co-linear with MQTLs on chromosomes 5H (MQTL5H.2), 2H (MQTL2H.3) and 3H (MQTL3H.1) in barley, respectively (Table 2). The ortho-MQTL for MQTL9-3, MQTL5-6, and MQTL9-3 on foxtail millet chromosomes was detected on chromosomes 4, 5, and 9 in barley (Table 2). Three foxtail millet MQTLs (MQTL1-5, MQTL9-1, MQTL9-3) were located in a syntenic region with detected MQTLs on rice chromosomes 2, 3, and 7, respectively, and foxtail millet MQTL2-1 was located in a syntenic region with detected MQTLs on rice chromosomes 1 and 7 (Table 2). Table S8 details all the ortho-MQTLs-localized genes and their respective annotations.
Table 2
Ortho-MQTLs in rice, barley and maize based on the syntenic region with identified MQTLs in foxtail millet
Ortho-MQTLs | Foxtail millet chr.no. (Genomic position in Mb) | Original MQTL name | Chr. no. (Genomic position in Mb) | References |
Barley |
ortho-MQTL2-5 | Chr. 2 (37.62–40.16) | MQTL5H.2 | Chr. 5H (527.14-571.19) | (Khahani et al. 2019) |
ortho-MQTL3-1 | Chr. 3 (1.01–1.92) | MQTL2H.3 | Chr. 2H (723.1-730.61) |
ortho-MQTL5-7 | Chr. 5 (45.22–45.62) | MQTL3H.1 | Chr. 3H (608.25-629.71) |
ortho-MQTL9-3 | Chr. 9 (6.19–7.01) | MQTL4H.2 | Chr. 4H (20.17–35.78) |
Maize |
ortho-MQTL5-6 | Chr. 5 (42.46–42.86) | MQTL-ASI-3 | Chr. 3(169.75-178.23) | (Almeida et al. 2013b) |
ortho-MQTL9-3 | Chr. 9 (6.19–7.01) | MQTL-GY-1b | Chr. 1(275.98-285.27) |
Rice |
ortho-MQTL1-5 | Chr. 1 (35.46–36.37) | MQTL2GW.8 MQTL2.YLD7 | Chr. 2 (28.41–29.7) Chr. 2 (29.3-31.84) | (Khahani et al. 2020) |
ortho-MQTL2-1 | Chr. 2 (1.47–2.42) | MQTL1.YLD1 MQTL7.GW23 MQTL7.PH20 MQTL7.HD15 | Chr. 1 (3.5–4.44) Chr. 7 (2.21–2.560) Chr. 7 (12.78–14.95) Chr. 7 (12.78–14.95) |
ortho-MQTL9-1 | Chr. 9 (3.47–4.12) | MQTL3.PH13 | Chr. 3 (32.92–33.03) |
ortho-MQTL9-3 | Chr. 9 (6.19–7.01) | MQTL3.YLD11 | Chr. 3 (27.82–29.55) |
Candidate Genes Located In The Mqtl Regions
The MQTL regions with CI of less than 1 Mb were explored looking for the drought stress up-regulated genes (according to RNA-seq and microarray data), leading to the identification of 101 candidate genes (Table S13). These potential candidate genes were categorized into several GO terms, for example, transcription regulation (e.g., SiHOX24, SiWRKY24, SiWRKY45, SiWRKY51, SiWRKY75, SiWRKY89, SiAP2/ERF, SiZAT12), signaling (e.g., SiPFK2, SiARF8, SiZAT12, SiSTK8, SiCSC1, SiWAK8, SiNCX, SiCDPK3, SiCKX5), transporter (e.g., SiPUB23), hormone metabolism (e.g., SiILL6 and SiACO), functional protein (SiHSP), glutathione metabolism (SiGGT1) and response to abiotic stress like drought and salinity (e.g., SiPHOS32) (Table S13).
Validation Of Potential Candidate Genes Via Qrt-pcr
The expression pattern of nine drought-regulated candidate genes was examined through qRT-PCR to validate the MQTLs results (Fig. 5). The expression profile of candidate genes was evaluated in the two drought-contrasting genotypes (Ilam1 as a drought-tolerant genotype and Kerman3 as a susceptible genotype). Eight genes were up-regulated (relative gene expression ranged between 6.19 and 1.74), and one gene was down-regulated (relative gene expression 4.28). A highly significant correlation was found between the MQTLs results and real-time RT-PCR profiling in Ilam1 genotype (R2 = 0.85) and Kerman3 genotype (R2 = 0.79).