Revealing the Genetic Architecture of Yield-Related and Quality Traits in Indian Mustard [ Brassica juncea (L.) Czern. and Coss.] Using Meta-QTL Analysis

: A meta-QTL analysis was conducted in Indian mustard to identify robust and stable meta-QTLs (MQTLs) by utilizing 1504 available QTLs, which included 891 QTLs for yield-related traits and 613 QTLs for quality traits. For yield-related traits, a total of 57 MQTLs (YRTs_MQTLs) were uncovered from the clustering of 560 projected QTLs, which had a 4.18-fold smaller confidence interval (CI) than that of the initial QTLs, whereas, for quality traits, as many as 51 MQTLs (Quality_MQTLs) were derived from 324 projected QTLs, which had a 2.65-fold smaller CI than that of the initial QTLs. Sixteen YRTs_MQTLs were observed to share chromosomal positions with 16 Quality_MQTLs. Moreover, four most promising YRTs_MQTLs and eight Quality-MQTLs were also selected and recommended for use in breeding programs. Four of these selected MQTLs were also validated with significant SNPs that were identified in previously published genome-wide association studies. Further, in silico functional analysis of some promising MQTLs allowed the detection of as many as 1435 genes, which also involved 15 high-confidence candidate genes (CGs) for yield-related traits and 46 high-confidence CGs for quality traits. After validation, the identified CGs can also be exploited to model the plant architecture and to improve quality traits through marker-assisted breeding, genetic engineering, and genome editing approaches


Introduction
Allotetraploid Indian mustard [Brassica juncea (L.) Czern and Coss] (AABB, 2n = 36) was synthesized via a number of hybridization events between B. rapa (AA, x = 10) and B. nigra (BB, x = 8) [1].In India, mustard is a widely grown oilseed crop, with an acreage of 7.2 million hectares and production of 7.1 million metric tonnes recorded during 2018-plant height, siliques per plant, seeds per siliqua, thousand-seed weight, number of primary and secondary branches, and siliques in the main shoot that contribute to increased seed yield, as well as the relationships between them, to select superior genotypes [2].
Seed oil, protein, and glucosinolates are the major factors that affect seed quality in B. juncea.The oil quality of Brassica species is determined by the relative proportions of the following unsaturated fatty acid fractions: 10-15% oleic acid (C18:1), 10-15% linoleic acid (C18:2), 14-16% linolenic acid (C18:3), and 50-60% erucic acid (C22:1) [3].Breeding efforts give preference to cultivars with high erucic acid content, which is needed for industrial purposes, and low erucic acid content, which is needed for human consumption [4].The increased oleic acid fraction in oil improves its appropriateness as a cooking medium [5], but the high level of linolenic acid is an undesirable feature that affects the oil's odor and taste, as well as its susceptibility to oxidation, limiting its shelf life [6].Breeding objectives aiming towards the improvement of these traits could be complicated owing to the pleiotropic associations among the traits, which have either negative or positive effects on other crucial yield-related traits [7,8].
Both yield-related and quality traits in B. juncea are inherited quantitatively, controlled by a large number of genes, and influenced by environmental factors.Over the last 25 years, more than 800 QTLs for yield-related traits and over 600 QTLs for quality traits have been identified (Table S1).Many factors, including the experimental conditions, size, and type of the mapping population, the density of genetic markers, statistical approaches, and the existence of genotype-environment and epistatic interactions, all affect the validity of QTL mapping findings [9][10][11].Moreover, the genetic intervals of these QTLs are often quite large, and the transfer of desired QTLs through marker-assisted breeding (MAB) may create the problem of linkage drag.The most important prerequisite for MAB and QTL cloning success is confining QTL to narrow genetic intervals and having a detailed understanding of the genetic regulation of a trait.Multiple populations tested at different locations and years may be used to achieve reliability in QTL work.Major consensus QTLs with consistency across populations and environments are considered most suitable for MAB [12,13].
Meta-QTL analyses have also been performed for some traits, including yield-related traits [12,37] and quality traits [13,38], in B. juncea by using QTL data derived from two to eight mapping populations phenotyped during the same study.These studies did not include QTL results from any of the previously published QTL mapping studies in their meta-analysis.Moreover, for these previously published meta-QTL studies on yieldrelated traits and quality traits in B. juncea, the maximum number of QTLs used was 187, which resulted in the detection of a maximum of 20 MQTLs for the studied traits [12].Several studies on QTL analysis for yield-related and quality traits in B. juncea have been conducted in the last two decades, and they have not been utilized for the prediction of MQTLs and detection of underlying CGs.
In this study, the QTL data associated with different yield-related and quality traits retrieved from 20 studies were used for the meta-QTL analysis.The main goals of this study were: (i) to create a highly dense consensus genetic map; (ii) to identify the most consistent and robust MQTLs with lower CIs than the initial QTLs; (iii) to provide accurate markers flanking the MQTLs for use in MAB; (iv) to identify the putative CGs within the promising MQTL regions; and (v) to conduct an expression analysis of the identified CGs.We think that our efforts will certainly help to increase the effectiveness of selection for different yield-related and quality traits in B. juncea.

QTL Data Collection for Meta-QTL Analysis
Using appropriate keywords, a comprehensive bibliographic search was conducted on Google Scholar (https://scholar.google.com/;accessed on 15 January 2022) and PubMed (http://www.ncbi.nlm.nih.gov/pubmed;accessed on 15 January 2022) for QTLs associated with different yield-related and quality traits reported during 1997-2021 (Table 1 and Table S1).The size of the mapping populations varied from 88-184 progenies of various types, which were mainly recombinant inbred lines (RILs) and doubled haploid (DH) populations phenotyped at different locations, each with more than one year.The following data were extracted from each study: (i) associated trait, (ii) linkage group/chromosome, (iii) genetic positions of the QTLs, (iv) size and type of bi-parental populations utilized, (v) logarithm of odds (LOD) score, (vi) flanking markers or closely linked markers, and (vii) PVE (phenotypic variation explained) or R 2 value.For the analysis, each QTL was given a unique identity based on the source studies, trait abbreviations, and linkage groups involved; QTLs for the same trait mapped on a single linkage group in an individual study were discriminated by numerical identifiers.
When data on the peak location were unavailable from a source paper, the most probable positions of the QTLs were calculated as the midpoint of the two markers associated with the QTLs.In some cases in which the LOD score was not specified, the likelihood ratio statistics given in the relevant source paper were used to estimate it (using the formula, LOD score = likelihood ratio statistics/4.6),and in others where it was simply stated that a LOD score of 3 was taken as a threshold to map the QTLs, the LOD of 3.0 was assumed for the concerned study.QTLs for different parameters were accommodated in the following two trait categories (Table S2): (i) yield-related traits, including siliquarelated traits (SRTs), seed number (SN), seed weight (SW), plant height (PH), number of branches (BN), and days to flowering (DTF), and (ii) quality traits, including protein content (PC), oil content (OC), erucic acid (C22:1), linolenic acid (C18:3), glucosinolate content (GSL), % pentyls, % butyls, % propyls, total aliphatics (TA), palmitic acid (C16:0), oleic acid (C18:1), eicosenoic acid (EA), linoleic acid (C18:2), stearic acid (C18:0), sinigrin (SIN), and gluconapin (GNA).
Some shared markers were also present in these genetic maps (reported in different studies), albeit at different genetic positions (Figure S1); for instance, for linkage group A1, eight markers were common between the maps derived from two different mapping populations reported in Jagannath et al. [44], and 26 markers were common between the maps retrieved from Lionneton et al. [46] and Lionneton et al. [48].Similarly, for A2, six markers were common between the maps derived from two different mapping populations reported in Jagannath et al. [44], and 15 markers were common between the maps reported in Lionneton et al. [46] and Lionneton et al. [48].Similarly, for A3, genetic maps from Dhaka et al. [37] and Ramchiary et al. [39] shared two markers, maps from Ramchiary et al. [39] and Yadava et al. [12] shared only one marker, maps from Lionneton et al. [46] and Lionneton et al. [48] shared 15 markers, and maps from Rout et al. [13] and Rout et al. [38] shared only one marker.
Similarly, for A4, one marker was found to be shared by maps derived from Yadava et al. [12] and Ramchiary et al. [39], and for A6, two markers were found to be shared by maps derived from Yadava et al. [12] and Ramchiary et al. [39].Similarly, maps obtained from Yadava et al. [12] and Ramchiary et al. [39] shared two markers, maps obtained from Yadava et al. [12] and Jagannath et al. [44] shared three markers, maps derived from two different mapping populations reported in Jagannath et al. [44] shared five markers, and maps obtained from Rout et al. [13] and Rout et al. [38] shared only one marker.Similarly, for A8, maps from Yadava et al. [12] and Ramchiary et al. [39] shared two markers, maps from Yadava et al. [12] and Jagannath et al. [44] shared one marker, maps from Jagannath et al. [44] and Khattak et al. [45] shared two markers, and maps from Rout et al. [13] and Rout et al. [38] shared two markers.Similarly, for A9, one marker was found to be shared by maps reported by Yadava et al. [12] and Ramchiary et al. [39], and six markers were found to be shared by maps derived from two different mapping populations reported by Jagannath et al. [44].
Similarly, for A10, maps from Dhaka et al. [37] and Aakansha et al. [40] shared one marker, maps from Yadava et al. [12] and Ramchiary et al. [39] shared one marker, and maps from Yadava et al. [12] and Jagannath et al. [44] also shared one marker.Similarly, only one marker was found to be shared by the Rout et al. [13] and Rout et al. [38] maps for B3.Similarly, four markers were reported to be shared between maps reported by Rout et al. [13] and Rout et al. [38] for B6.Similarly, for B8, only one marker was found to be shared by the maps reported by Rout et al. [13] and Aakansha et al. [40].However, no shared markers were found between individual A5, B1, B2, B4, B5, and B7 maps retrieved from different studies.
The R package LPmerge (Endelman, Madison, WI, USA) was utilized for the development of the consensus map needed for the meta-analysis [54].LPmerge creates a consensus map by integrating genetic maps from different populations.This software utilizes linear programming to effectively reduce the mean absolute error between the genetic maps from individual populations and the consensus map.This minimization is carried out under the constraints of linear inequality, ensuring that the order of the marker loci in the genetic maps is preserved.When linkage maps have incompatible marker orders, a minimum set of ordinal restrictions is removed to resolve the problems [54].Overall, there were two steps in the process of creating a consensus map using LPmerge: The first step provided statistics on the numbers of consensus bins, markers, and ordinal conflicts, and the second step provided the root-mean-square error (RMSE) values between each initial map and the consensus map, along with the mean and standard deviation of the RMSE.This program generates consensus maps for K = 1 to 4, where K is the largest possible interval size.For the meta-analysis, a map with a length that was similar to the mean of the component genetic maps and the lowest mean and standard deviation for RMSE was chosen.

QTL Projection
Details on QTLs (viz., genetic positions, R 2 /PVE values, and LOD scores) associated with the studied traits were compiled and projected separately for the yield-related and quality traits groups onto the newly generated consensus map.QTL data with no data on genetic positions, PVE values, etc. were not utilized for the projection.However, in some cases, when peak positions were available but CIs were not given, the CI (95%) was calculated using the following population-specific equations: where, N is the number of the individuals of the population utilized for mapping and PVE is the contribution of an individual QTL to the total trait variation.The values 287 and 163 are the constants derived from the different simulations [22,55].Equation ( 1) was used for the QTLs mapped using DH populations [55] and Equation ( 2) was used for the QTLs mapped using RIL populations [22].The QTLProj tool furnished in the software BioMercator v4.2 [9,56,57] was used for QTL projection.QTLProj uses a dynamic method to choose the ideal projection context.A pair of common markers that flank the QTL in the original map (provided in individual mapping studies) and have consistent distance estimates between the maps make up an ideal context.Following two parameters determine how the algorithm behaves when trying to find such a configuration: (i) the minimal ratio of the flanking marker interval distances and (ii) the minimal p-value discovered when examining the homogeneity (or similarity) of the flanking marker interval distances between the initial and consensus maps [56,57].

Prediction of MQTLs Associated with Yield-Related and Quality Traits and Their Validation with Genome-Wide Association Studies
The meta-QTL analysis was carried out using the BioMercator v4.2 software, and it was performed separately for each group of traits with two different methods: (a) Veyrieras's two-step algorithm [58] and (b) Goffinet and Gerber's approach [10].The Veyrieras two-step algorithm was used for the meta-QTL analysis when more than 10 projected QTLs were available on a single linkage group.The best meta-QTL model was chosen based on the most common value among the Akaike Information Criterion (AIC), corrected AIC, Bayesian Information Criterion (BIC), and Average Weight of Evidence (AWE) criteria.However, Goffinet and Gerber's method was used when the total number of QTLs projected on a given linkage group was less than 10.Following that, the consensus QTLs discovered using the best model (that with the lowest AIC value) were considered MQTLs.The mathematical methods and algorithms used in BioMercator v4.2.3 (Yannick De Oliveira, Paris, France) have been well documented [9,56,57].In addition to the abovementioned criteria, only genomic regions with at least two QTLs (each from a different study) co-localized were considered MQTLs.Furthermore, a search was conducted for previously published genome-wide association (GWA) studies on yield-related and quality traits in Indian mustard.The data on significant SNPs and their genomic positions were collected from the collected studies.The genomic positions of identified significant SNPs were compared to the physical positions of MQTLs detected in individual linkage groups.MQTLs that co-localized with at least one significant SNP were considered GWAS validated/verified.

Detection of Candidate Genes and Their Expression Analysis
Some most stable and robust MQTLs were also considered for the detection of putative candidate genes (CGs).The first step was to extract the nucleotide sequences of the flanking markers of chosen MQTLs from the Brassica genome database (https://www.brassica.info/;accessed on 15 May 2022) or the available literature.BLAST searches were then performed using the retrieved nucleotide sequences against the reference genome of B. juncea (Braju tum V 2.0) available on the Brassica database (http://brassicadb.cn/#/;accessed on 20 May 2022) to get the genomic regions of the corresponding MQTLs.MQTLs occupying genomic regions of less than 1 Mb were directly investigated for the detection of CGs using the Brassica database.However, for MQTLs with physical intervals of more than 1 Mb, first, the peak physical positions of MQTLs were calculated by using a formula recently proposed by us [29], which is as follows: Peak position(bp) = start position(bp) + ((end position(bp Then, the total 1 Mb genomic region of each such MQTL (0.5 Mb flanking region of the MQTL peak) was explored for the detection of available CGs. Predicted genes were then annotated by BLAST analysis against the Arabidopsis thaliana database (http://brassicadb.cn/#/Annotations/;accessed on 23 May 2022).Further, different developmental gene expression datasets were downloaded from BioAnalytic Resource in the e-Northerns w.Expression Browser (http://bar.utoronto.ca;accessed on 12 June 2022), and the array was normalized with a TGT value of 100 using the Gene Chip Operating Software.For expression analysis, datasets comprised of different samples of relevant developmental stages (details are provided in Schmid et al. [59]) were downloaded, and then the differential expressions of some selected genes were shown through a heatmap using the default criteria of the webtool 'Clustvis' (https://biit.cs.ut.ee/clustvis/; accessed on 16 June 2022; [60]).

Distribution of QTLs on the B. juncea Genome and Their Salient Features
A total of 1504 QTLs, including 891 QTLs for yield-related traits and 613 QTLs for quality traits assessed under non-stress conditions, were collected from twenty mapping studies (Table S1). Figure 1 depicts salient features of the QTLs associated with different yield-related traits, such as the number of QTLs per trait (Figure 1a), number of QTLs per chromosome (Figure 1b), PVE of individual QTLs (Figure 1c), LOD scores of individual QTLs (Figure 1d), and CIs of the initial QTLs (Figure 1e).Of the total QTLs that were collected for yield-related traits, 57.7% (514/891) of the QTLs were found in sub-genome A, while 42.3% (377/891) were found in sub-genome B.More than 98% of the total QTLs (880/891) had data needed for meta-analysis.The PVE (%) by a single QTL varied between 2.5 and 50.1%, with an average of 9.45%.The CIs of QTLs ranged from 0 to 37 cM, with a mean of 12.40 cM.With a total of 113 QTLs, chromosome A7 carried the highest number of QTLs for the studied parameters, followed by chromosome B3 (103 QTLs) and chromosome A3 (80 QTLs), while chromosome B6 carried the minimum (12 QTLs).Similarly, Figure 2 presents major characteristics of the QTLs associated with quality traits, such as the number of QTLs per trait (Figure 2a), chromosome-wise distribution of QTLs (Figure 2b), PVE of individual QTLs (Figure 2c), LOD scores of individual QTLs (Figure 2d), and CIs of the initial QTLs (Figure 2e).
Of the total QTLs that were collected for quality traits, 57.8% (354/613) of the QTLs were carried by sub-genome A, and 42.2% (259/613) QTLs were carried by sub-genome B.More than 80% (499/613) of the collected QTLs had data needed for meta-analysis.The PVE (%) by an individual QTL varied from 1 to 92%, with an average of 16.83%.The CIs of QTLs ranged from 0 to 53 cM, with an average of 11.44 cM.With 75 QTLs, chromosome A2 had the maximum QTLs for the parameters studied, followed by chromosomes A3 (67 QTLs) and A8 (60 QTLs), whereas chromosomes A4 and A10 had the minimum (each with 11 QTLs).

Features of the Consensus Map
The consensus map involved a total of 13,725 markers of different types, including RFLP, AFLP, SSR, CAPS, SNPs, IP, RAPD, and some gene-specific markers (Table S3).The total length of the consensus map measured 3493.8 cM, with an average chromosomal length of 194.1 cM and individual chromosomal lengths ranging from 132 (A6) to 299.5 cM (A5).The average number of markers per chromosome was around 762, with the highest number of markers (1343 markers) on chromosome B3 and the lowest number (308) on chromosome B4 (Table S4).The whole-genome marker density was 3.93 markers/cM, with a minimum (1.73 marker/cM) on chromosome B4 and a maximum (6.39 markers/cM) on chromosome B1.Sub-genome A had 7394 markers over a distance of 2019.8 cM, while sub-genome B had 6331 markers over a distance of 1474 cM.

MQTLs for Yield-Related Traits and Their Distribution in the Genome
Seven hundred and ninety-eight QTLs associated with yield-related traits were successfully projected onto the consensus map, resulting in 57 MQTLs (each carrying at least two QTLs from different studies) (Table S5, Figure 3a) derived from 560 QTLs, while nine remained as singletons (i.e., single QTLs), and three QTLs were outside the MQTL supporting intervals.Of the 798 QTLs, 226 QTLs were clustered into 38 genomic regions, which were termed as QTL hotspots (each carrying multiple QTLs from only a single study) (Figure 3b).Ninety-three QTLs could not be projected owing to either of the reasons discussed in our earlier reports [33][34][35].Thirty-one MQTLs were found in subgenome A and 26 were found in sub-genome B. The lowest number of QTLs (10 QTLs) were projected on chromosome B6, while the maximum number of QTLs (99 QTLs) were projected on chromosome A7, which yielded five MQTLs.Overall, the number of MQTLs ranged from one on chromosome B5 to a maximum of six MQTLs each on chromosomes A3 and B3, with no MQTLs identified on chromosomes A1 and A5.Twenty MQTLs were obtained from the grouping of at least 10 initial QTLs involving different mapping populations.The number of traits/parameters associated with individual MQTLs varied from a minimum of one (YRTs_MQTLA2.2) to 12 (YRTs_MQTLA10.1)(Table S5).The LOD value of the MQTLs associated with yield-related traits varied from 2.78 to 16.28 (with a mean of 4.78), whereas the PVE value varied from 5.38 to 24.63% (with a mean of 8.82%) (Table S6).There were 42 YRTs_MQTLs that involved at least one major QTL (with PVE ≥10%) (Table S7), whereas the remaining 15 YRTs_MQTLs included only the minoreffect QTLs (with PVE <10%).Further, 21 QTL hotspots were found in sub-genome A and 17 were found in subgenome B. The lowest number of QTLs (only two QTLs) were projected on chromosome A7, while the maximum number of QTLs (43 QTLs) were projected on chromosome A1, which gave the origin to five QTL hotspots.Overall, the number of QTL hotspots ranged from one each on chromosomes A4, B2, and B7 to a maximum of five hotspots each on chromosomes A1, A5, and B8, with no hotspots identified on chromosomes A8, A10, B1, and B6.The LOD values of the QTL hotspots associated with yield-related traits ranged from 2.91 to 10.5 (with a mean of 4.48), whereas the PVE values ranged from 4.69 to 11.45% (with a mean of 7.52%).
The CIs of the MQTLs and QTL hotspots varied from 0.15 to 9.05 cM, with an average of 3.03 cM, demonstrating a substantial reduction from the projected QTLs, which had CIs ranging from 0 to 37 cM, with a mean of 12.66 cM (Figure 3c).In addition, each of the 18 MQTLs and the 13 QTL hotspots had a CI of ≤2 cM.Overall, there were considerable differences between different chromosomes, with the average CIs of MQTLs and QTL hotspots being 4.18-fold lower than those of the projected initial QTLs (Figure 3c).On chromosomes A3 and B2, the average CIs of MQTLs and QTL hotspots decreased by 7.15 and 6.92 times, respectively.On chromosomes A7 and A10, the average CIs decreased by 5.43 and 5.24 times (Figure 3c).In addition, Table S7 shows the CIs of the original major QTLs associated with yield-related traits, as well as the fold reduction in the CIs of individual major QTLs after the meta-analysis.The number of traits/parameters associated with QTL hotspots varied from a minimum of one (YRTs_HotspotA2.2) to nine (YRTs_HotspotA1.5).Table S8 includes detailed information on the QTL hotspots.

MQTLs for Quality Traits and Their Distribution in the Genome
Out of the 613 total QTLs associated with quality traits, 410 were successfully projected onto the consensus map, resulting in 51 MQTLs derived from 324 QTLs, while 18 remained as singletons and two QTLs were outside the MQTL supporting intervals (Table S9, Figure 4a).There were also 66 QTLs among the 410 projected QTLs, which were grouped into 20 QTL hotspots.For either of the reasons mentioned in our earlier reports of meta-analysis [33][34][35], two hundred and three QTLs were not projected.There were 28 MQTLs in sub-genome A and 23 in sub-genome B. On chromosome B5, the lowest number of QTLs (five QTLs) were projected, whereas the largest number of QTLs (47 QTLs) were projected on chromosome B7, which gave rise to four MQTLs.Overall, there were two MQTLs on each of chromosomes A1, A6, and B8, with a maximum of six MQTLs on chromosomes A2, B1, and B2, and no MQTLs on chromosomes A4, B4, and B5.Twelve MQTLs were obtained from the grouping of at least 10 initial QTLs involving different mapping populations (Table S9).The number of traits/parameters associated with individual MQTLs varied from a minimum of one in eight MQTLs (Quality_MQTLA8.3,A10.1, A10.2, A10.3, B1.2, B2.4,B3.2, B7.2) to nine in one MQTL (Quality_MQTLA9.2).The LOD values of the individual MQTLs associated with quality traits ranged from 2.65 to 39.79 (with an average of 7.6), whereas the PVE values ranged from 4.12 to 33.08% (with a mean of 12.52%) (Table S10).Among the 51 Quality_MQTLs, 39 Quality_MQTLs involved at least one major QTL (with PVE ≥10%) (Table S7), and the remaining 12 Quality_MQTLs included only the minor-effect QTLs.
In addition, 14 QTL hotspots in sub-genome A and six in sub-genome B were discovered (Figure 4b).The lowest number of QTLs (only two) were projected on chromosomes A2 and A5, whereas the maximum number of QTLs (15) were projected on chromosome A3, which gave rise to four QTL hotspots.Overall, there was one QTL hotspot on each of chromosomes A2, A5, A7, B4, and B7, and a maximum of four hotspots on chromosome A3, with no hotspots found on chromosomes A1, A6, A10, B1, B2, B5, and B8.Individual QTL hotspots associated with quality traits had LOD values ranging from 3.25 to 35.87 (with a mean of 6.52) and PVE values ranging from 2.7 to 45.67% (with a mean of 11.39%).The CIs of the MQTLs and QTL hotspots ranged from 0.19 to 33.41 cM, with an average of 4.43 cM, demonstrating a substantial reduction from the projected QTLs, which had CIs ranging from 0 to 53 cM, with a mean of 11.76 cM.Furthermore, each of the 18 MQTLs and six QTL hotspots had a CI of less than 2 cM (Figure 4c).Overall, the average CI of MQTLs and QTL hotspots was 2.65 times lower than that of the projected QTLs, with significant differences among the chromosomes.Table S11 includes information on the CIs of the original major QTLs associated with quality traits, as well as the fold reduction in the CIs of individual major QTLs following the meta-analysis.The number of traits/parameters associated with individual QTL hotspots varied from a minimum of one in nine hotspots to a maximum of only two traits in the remaining 11 hotspots.Table S12 includes detailed information on the QTL hotspots associated with quality traits.The distribution of all YRTs_MQTLs and Quality_MQTLs is depicted in Figure 5.
A search for putative genes of the above five YRTs_MQTLs and eight Quality_MQTLs yielded a total of 1435 genes (including 436 genes available from YRTs_MQTLs and 999 from Quality_MQTLs) in B. juncea (Table S14) with the corresponding orthologs in A. thaliana genome (Table S15).These genes included 232 genes (52 genes from YRTs_MQTLs and 180 genes from Quality_MQTLs) with no function descriptions available.Among the YRTs_MQTLs, the maximum number of 344 genes were associated with YRTs_MQTLA7.2,and a minimum of seven genes were available in YRTs_MQTLA10.1.However, among the Quality_MQTLs, the maximum number of 185 genes were associated with Quality_MQTLA8.2 and a minimum of 106 genes were available in Quality_MQTLA8.1.Several genes/gene families with similar proteins/products were repeatedly detected in YRTs_MQTL regions, for instance, 18 genes for protein kinases, eight for glycoside hydrolases and zinc finger proteins each, six for glycosyl transferase and SANT/Myb-domain-containing proteins, etc.Similarly, several genes with similar proteins were repeatedly identified in Quality_MQTL regions, for instance, 38 genes for zinc finger containing proteins, 26 for protein kinases, 15 for glycoside hydrolases, 13 for SANT/Myb-domain-containing proteins, 11 for AP2/ERF and F-box-domain-containing proteins each, and eight for glycosyl transferases and oxoglutarate/iron-dependent dioxygenases each.
Further, based on the available literature and functional annotations of A. thaliana orthologs, 15 high-confidence CGs (two CGs corresponded to a single ortholog in A. thaliana, leaving 14 unique CGs) associated with yield-related traits and 46 highconfidence CGs (two CGs corresponded to a single ortholog in A. thaliana, leaving 45 unique CGs) associated with quality traits were selected (Table 2).Candidate genes from YRTs_MQTLs included those encoding for various transcription factors (TFs), including AP2/ERF, MADS box, SANT/Myb, Tify, and zinc-finger-domain-containing proteins, protein kinases, glycoside hydrolases, glutathione S-transferases, etc., whereas Quality_MQTLs included those encoding for the following proteins/products: acyl-CoA oxidases, enoyl-acyl-CoA oxidases, enoyl-CoA hydratase/isomerases, FAD-linked oxidases, fatty acid hydroxylase, fatty acyl-coenzyme A reductases, GDSL lipase/esterases, oleosins, oxoglutarate/iron-dependent dioxygenases, palmitoyltransferases, UDP-glucuronosyl/UDP-glucosyltransferases, very-long-chain 3ketoacyl-CoA synthases, different TFs (e.g., WRKY, MADS box, SANT/Myb), and protein kinases.An expression analysis of 14 orthologs of CGs available from YRTs_MQTLs (Figure 6) and 42 orthologs of CGs available from Quality_MQTLs (Figure 7) was conducted.Among the CGs associated with YRTs_MQTLs, some of the genes showed comparatively high expressions at the vegetative stage (e.g., leaves, shoots, vegetative rosettes at different time intervals); for instance, the genes encoding for the following products/proteins: Tify-TF, protein kinase, zinc-finger-domain-containing protein, glycoside hydrolase, etc.However, some genes showed high expressions mainly at the reproductive stage (e.g., young and developed flowers, shoot apex, siliques, seeds); for instance, the genes encoding for the following products/proteins: major intrinsic protein, glutathione S-transferase, SANT/Myb-domain-containing proteins, etc.Some of the genes, such as those encoding for the homeobox domain, AP2/ERF TF, MADS-box TF, protein kinase, etc., showed almost consistent expressions across all the developmental stages studied (Figure 6).Similarly, among the CGs associated with Quality_MQTLs, some CGs showed consistent expressions in almost all the siliques and seed-related tissues (Figure 7); for instance, glycoside hydrolase, fatty acid hydroxylase, oleosin, very-longchain 3-ketoacyl-CoA synthase, acetyl-coenzyme A carboxyltransferase, and NADH dehydrogenase.

Discussion
Breeding for higher yield and seed quality traits is crucial for increasing the productivity of the Indian mustard crop.Meta-QTL analysis effectively eliminates the influence of population type and size, genetic background, and phenotyping environment on QTLs, and it integrates QTL data derived from multiple independent studies [9,12].This method has been widely applied in many crop plants, such as wheat [29,30,[33][34][35], rice [14,15], maize [18,19], and soybean [22][23][24][25][26], including those in B. napus [61,62] and B. juncea [12,19,37,38], for numerous traits.This study could potentially be used in Indian mustard to identify QTL clusters in defined regions in a large set of yield-related and quality trait data obtained from twenty populations.

MQTLs for Yield-Related and Quality Traits
During the present study, we used hundreds of initial QTLs, thereby ensuring more systematic and reliable anchoring of genetic loci, as the consistency of statistical findings has been found to be significantly and positively associated with the number of initial QTLs utilized for meta-analysis [29].These QTLs were associated with different traits, which were grouped into two major trait categories-(i) yield-related traits and (ii) quality traits-for better representation and interpretation of the results (individual traits associated with each MQTL are also provided; see Tables S5 and S9).The collected QTLs were unevenly distributed on the chromosomes, consistently with previous studies; subgenome A carried more QTLs for both quality (354 QTLs) [13] and yield-related traits (514 QTLs) [12].Of the 798 total QTLs collected, 560 QTLs were utilized for the identification of 57 MQTLs for yield-related traits, whereas, for quality traits, 324 out of 613 QTLs were utilized for the detection of 51 MQTLs, implying that the current study is, so far, the most comprehensive meta-QTL analysis for the studied traits in B. juncea.
The distribution of MQTLs and initial QTLs on different chromosomes was not consistent, which was primarily owing to the different numbers of initial QTLs clustered in MQTLs.In this study, 35.08% (20/57) of the MQTLs associated with yield-related traits were composed of a minimum of 10 initial QTLs, and one (YRTs_MQTLA7.2) of them was composed of as many as 53 initial QTLs.As many as 12 (24% of the total) MQTLs associated with quality traits were composed of a minimum of 10 initial QTLs, and four of them were composed of at least 15 QTLs, which was significantly higher than that found in previous studies [12,13,37].These predicted MQTLs that contained initial QTLs identified from different bi-parental populations are considered more accurate and stable MQTLs.Following the criteria laid down by Löffler et al. [63] and later improved by Saini et al. [29], we selected at least four 'breeder's MQTLs' for yield-related traits, viz., YRTs_MQTLA7.2,A7.3, A9.2, and B3.1 (possessed <2 CI, 10 to 53 initial QTLs were involved in the development of these MQTLs, LOD score for a single MQTL varied from 5.58 to 9.4, and PVE value for a single MQTL varied from 10.40 to 17.24%) (Table 3) and eight breeder's MQTLs for quality traits, which included the following: Quality_MQTLA2.4,A2.5, A3.1, A3.2, A8.1, A8.2, B7.3, and B7.4 (possessed <2 CI, 10 to 29 initial QTLs were involved in the formation of these MQTLs, LOD score for a single MQTL ranged from 12.86 to 39.79, and PVE value for a single MQTL ranged from 13.64 to 33.08%) (Table 4).
In addition to QTL mapping, a genome-wide association study (GWAS) based on high-throughput markers is another high-precision method for identifying genomic regions associated with complex quantitative traits in crop plants [64].Further, the combination of GWAS and meta-QTL analyses can effectively integrate the significant SNPs (identified through GWAS) and QTLs (identified through linkage-based interval mappings) from different studies and help better understand the genetic architecture of complex quantitative traits [27,29,32,36,65,66].However, compared to other major crops, only a few GWA studies are available for Indian mustard, and these include three of our own studies, those of Akhtar et al. [67][68][69], and three other studies, each of which was published in a different country, namely, India [70], Australia [71], and China [72].Further, two of these GWA studies [67,72] did not provide the complete information (such as physical positions of significant SNPs) required for the validation of the MQTLs.
These selected breeder's MQTLs could also be validated before being utilized in breeding programs.Traditional QTL validation requires the development of near-isogenic lines to assess the effects of detected QTLs on the target trait.Because of its inclusive QTL property, MQTL validation differs intuitively from traditional QTL validation.The identification of MQTL in a condensed genomic region according to the algorithmic definition of the meta-analysis ensures the validation of known QTLs in the region.As a result, identifying and validating a linked surrogate marker on a set of genotypes with contrasting phenotypes may be practically important.The presence of at least a few QTL alleles in the MQTL region may be ensured by the significant association of marker alleles with phenotypes.For YRTs_MQTLs, the linked markers can be validated using two different sets of genotypes that are contrasted for the target yield component traits or yield per se, whereas for quality MQTLs, the linked markers can be validated using different sets of genotypes that are contrasted for the target quality trait(s) or overall quality.Further, the QTL hotspots identified for different yield-related and quality traits in the present study may be considered as population-specific QTLs or unique QTLs, which may also be investigated further for the traits under study.5), SD (7), PB (8), DTF (10), PH (10), MSL (7), SB a For abbreviations, see Table 1.
a For abbreviations, see Table 1.

Co-Localization of MQTLs: Inferring a Genetic Relationship between Yield-Related and Quality Traits
Genetic improvements in seed yield and quality traits, including oil content, fatty acid composition, and protein content, are the central targets in mustard breeding programs.Analyses of the genetic relationships of quantitative traits have been important for evaluating the applicability of joint selection of two or more traits [73].In the present study, sixteen MQTLs associated with yield-related traits were found to have genetic positions that almost overlapped those occupied by 16 MQTLs associated with quality traits (Table S16, Figure 5).The relationship between two or more desired traits makes it easier for breeders to improve them at the same time.Co-localization of QTLs for yieldrelated and quality traits is in accordance with the findings of previous studies that revealed a correlation among these traits in B. juncea [7,12,13,46].Therefore, the markers flanking these overlapped MQTLs can be effectively utilized for simultaneous improvements of different sets of yield-related and quality traits through MAB.The presence of a tight linkage between genes in the same genomic regions controlling different traits or the pleiotropy of underlying genes may be another reason for this colocalization of MQTLs associated with different traits [13].It is also evident from genetic studies that the same genomic region may affect different traits in an opposite manner [7,74].The substantial negative relationships between different traits (e.g., seed oil and protein content), their common/linked QTLs, and epistasis with opposite effects provide potential problems for plant breeders attempting to improve these traits in Brassicas at the same time [74].

Candidate Genes for Yield-Related and Quality Traits
MQTLs are regarded as potential genomic regions that may contain the most promising CGs for the traits under study.A total of 1435 genes were found in the current study after gene mining was carried out in the five most robust and stable YRTs_MQTL regions and the eight Quality_MQTL regions.An important phenomenon is that several genes (e.g., protein kinases, UDP-glucuronosyl/UDP-glucosyltransferase, proteinase inhibitor I3, and glycoside hydrolase) appeared as gene clusters in the MQTL regions.Plant genomes contain many of these gene clusters, and the proteins they produce are involved in a wide range of enzymatic processes [29,75].Additionally, it was shown that members of a gene cluster might be found close together (a few thousand base pairs apart) in a small genomic region, thus encoding similar products/proteins and, therefore, sharing a generalized function.Furthermore, to investigate potential roles in the regulation of the traits in question, additional research should be carried out on the genes encoding hypothetical proteins and those lacking function descriptions.A total of 61 highconfidence CGs associated with the traits (15 with yield-related traits and 46 with quality traits) under study were also selected based on a survey of the available literature and functional annotation of corresponding A. thaliana orthologs.Owing to the lack of transcriptomic data repositories for B. juncea, we decided to use transcriptomic expression data for A. thaliana to explore the expression patterns of these high-confidence CGs at different developmental stages.Genes mainly expressed at vegetative stages may have their major roles in the regulation of yield-related traits, and genes expressed at reproductive stages may have their association with both yield-related and quality traits.However, the genes that showed almost consistent expressions across all of the developmental stages may have a pleiotropic effect on both yield-related and quality traits.
The functions of some of the high-confidence CGs associated with yield-related traits can be discussed as follows: (i) The TFs of the AP2/ERF family serve as crucial regulators in many physiological and biological processes, including plant morphogenesis, hormone signal transduction, response mechanisms to different stresses, and metabolite regulation [76].(ii) The ortholog of the BjuVA07G37930 gene, i.e., AGL12 encoding for MADS TF, in A. thaliana has been shown to have an important role in root development, as well as in flowering transition [77].(iii) The ortholog of the BjuVA07G37280 gene, i.e., JAZ9 encoding for Tify-domain-containing protein, regulates jasmonte hormone signaling, which plays an important role in controlling the growth and development of Arabidopsis plants [78].(iv) The ortholog of the BjuVA09G08320 gene, i.e., SHA1, is known to play a key role in shoot apical meristem maintenance in A. Thaliana [79].In higher plants, the primary shoot apical meristem is established during embryogenesis and contributes irregularly throughout the plant's lifetime to the production of all aerial organs, including leaves, stems, and flowers [79].(v) Earlier, a mitogen-activated protein kinase (i.e., OsMAPK6) was reported to affect grain size and biomass production in rice via its involvement in cell proliferation, brassinosteroid signaling, and homeostasis [80].Similarly, the following is a brief summary of the functions of some of the high-confidence CGs associated with quality traits: (i) Over-expression of Gly-Asp-Ser-Leu (GDSL)-motif esterase gene (i.e., EgGDSL) significantly improved total fatty acid content, drastically reduced the levels of linoleic acid and 11-eicosenoic acid, and increased the content of stearic acid in the seeds of transgenic lines [81].Functions of different members of the GDSL lipase family in improving the oil content have also been reported in some other studies [82].(ii) Overexpression of a gene encoding for WRKY TF (WRKY43) increased polyunsaturated fatty acid content in seeds of A. thaliana [83].(iii) The accumulation of oleosins regulates the size of seed oil bodies in A. thaliana [84].(iv) Recently, in 2017, Seo and Kim [85] discussed various roles of different myeloblastosis (Myb)-like TFs in glucosinolate biosynthesis in Brassicaceae.(v) FAD-linked oxidases are known to regulate the content of dienoic fatty acids and increase the resistance toward different abiotic stresses in plants [86].(vi) Verylong-chain 3-ketoacyl-CoA synthases are known to be involved in the biosynthesis of very-long-chain fatty acids in plants [87].Once validated, these CGs may be used as promising targets that can be exploited by using biotechnological approaches to improve the yield-related and quality traits in Indian mustard.

Conclusions
The meta-analysis of a large number of QTLs derived from independent studies provides useful knowledge for consensus genomic regions, i.e., QTLs, to be exploited for molecular breeding.During the present study, data were collected on as many as 1504 QTLs, including 891 QTLs for yield-related traits and 613 QTLs for quality traits, from studies published from 1997 to 2021.Using these QTLs, we identified 57 and 51 MQTLs for yield-related and quality traits, respectively.Among these, four very promising MQTLs (YRTs_MQTLA7.2,A7.3, A9.2, and B3.1) were recommended for the improvement of the yield of plants, whereas at least eight MQTLs (Quality_MQTLA2.4,A2.5, A3.1, A3.2, A8.1, A8.2, B7.3, and B7.4) were selected and recommended for use in breeding programs to improve quality traits.Further, in silico investigation of some selected promising MQTLs allowed the identification of as many as 1435 genes, including 15 high-confidence CGs associated with yield-related traits and 46 CGs associated with quality traits.Nevertheless, after any in silico analysis of this kind, in vivo confirmation/validation of any of the detected loci, specifically the CGs, is required, which may be achieved using any of the available approaches, such as gene cloning, gene silencing, and even other advanced approaches, such as transcriptomics, proteomics, and metabolomics.The results of this study provide a basis for understanding the genetic underpinnings of yield-related and quality traits in Indian mustard.In general, the MQTL and CG information presented in this study should be helpful for the establishment of functional markers and the introduction of elite alleles in mustard breeding efforts.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/agronomy12102442/s1. Table S1.Initial QTLs associated with yield-related and quality traits used in this study; Table S2.Trait categories considered in the present study; Table S3.Consensus map constructed in the present study; Table S4.Details of the consensus map constructed in the present study; Table S5.MQTLs associated with yield-related traits identified in the present study; Table S6.LOD scores and PVE values of individual QTLs involved in different YRTs_MQTLs; Table S7.Fold reduction in confidence intervals of major QTLs associated with yield-related traits after meta-analysis; Table S8.QTL hotspots associated with yield-related traits identified during the present study; Table S9.MQTLs associated with quality traits identified in the present study; Table S10.LOD scores and PVE values of individual QTLs involved in different Quality_MQTLs; Table S11.Fold reduction in confidence intervals of major QTLs associated with quality traits after meta-analysis; Table S12.QTL hotspots associated with quality traits identified during the present study; Table S13.GWA studies on yield-related and quality traits (used for comparing MQTLs identified in the present study); Table S14.Putative CGs identified within the promising MQTL regions; Table S15.Arabidopsis orthologs of the genes identified in the MQTL regions; Table S16.Co-localization of Quality_MQTLs with YRTs_MQTLs reported in the present study; Figure S1.Shared markers among the linkage maps derived from different interval mapping studies.

Figure 1 .
Figure 1.Salient features of the QTLs associated with yield-related traits: (a) number of QTLs per trait, (b) number of QTLs per chromosome, (c) phenotypic variation explained (PVE) by individual QTLs, (d) LOD scores of individual QTLs, and (e) confidence interval (CI) of the QTLs.Trait abbreviations: DTF: days to flowering, SRTs: siliqua-related traits, SN: seed number, SW: seed weight, PH: plant height, and BN: number of branches.

Figure 3 .
Figure 3. Salient features of YRTs_MQTLs and QTL hotspots: (a) number of QTLs involved in different MQTLs on different chromosomes, (b) number of QTLs involved in QTL hotspots mapped on different chromosomes, and (c) comparison of CIs of projected QTLs and CIs of MQTLs and QTL hotspots.

Figure 4 .
Figure 4. Salient features of Quality_MQTLs and QTL hotspots: (a) number of QTLs involved in different MQTLs on different chromosomes, (b) number of QTLs involved in QTL hotspots mapped on different chromosomes, and (c) comparison of CIs of projected QTLs and CIs of MQTLs and QTL hotspots.

Figure 5 .
Figure 5. Distribution of YRT_MQTLs (shown in red), Quality_MQTLs (shown in black), and selected breeder's MQTLs (shown in yellow and green) on different chromosomes.

Figure 6 .
Figure 6.Expression patterns of high-confidence CGs associated with yield-related traits in different plant tissues.

Figure 7 .
Figure 7. Expression patterns of high-confidence CGs associated with quality traits in different plant tissues.

Table 1 .
Summary of the mapping experiments utilized in the present study for meta-QTL analysis.

Table 2 .
High-confidence candidate genes available from some promising MQTLs.
a For abbreviations, see Table1.

Table 3 .
Breeder's MQTLs for yield-related traits identified in this study.

Table 4 .
Breeder's MQTLs for quality traits identified in this study.