Integration of QTL Mapping and Gene Fishing Techniques to Dissect the Multi-Main Stem Trait in Rapeseed (Brassica napus L.)

Rapeseed is one of the most important oilseed crops in the world. Improving the production of rapeseed is beneficial to relieve the shortage of edible vegetable oil. As the organ of support and transport, the main stem of rapeseed controls the plant architecture, transports the water and nutrients, and determines the number of inflorescence. Increasing the number of main stems would be helpful for the yield improvement in Brassica napus (B. napus). This attractive multi-main stem (MMS) trait was observed in the KN DH population. We investigated not only the frequency of MMS traits but also dissected the genetic basis with QTL mapping analysis and Gene-Fishing technique. A total of 43 QTLs were identified for MMS based on high-density linkage map, which explained 2.95–14.9% of the phenotypic variation, among which two environmental stable QTLs (cqMMS.A3-2 and cqMMS.C3-5) were identified in winter and semi-winter environments. Epistatic interaction analysis indicated cqMMS.C3-5 was an important loci for MMS. According to the functional annotation, 159 candidate genes within QTL confidence intervals, corresponding to 148 Arabidopsis thaliana (A. thaliana) homologous genes, were identified, which regulated lateral bud development and tiller of stem, such as shoot meristemless (STM), WUSCHEL-regulated-related genes, cytokinin response factors (CRF5), cytokinin oxidase (CKX4), gibberellin-regulated (RDK1), auxin-regulated gene (ARL, IAR4), and auxin-mediated signaling gene (STV1). Based on Gene-Fishing analysis between the natural plants and the double-main stem (DMS) plant, 31 differentially expressed genes (DEGs) were also obtained, which were related to differentiation and formation of lateral buds, biotic stimulus, defense response, drought and salt-stress responses, as well as cold-response functional genes. In addition, by combining the candidate genes in QTL regions with the DEGs that were obtained by Gene-Fishing technique, six common candidate genes (RPT2A, HLR, CRK, LRR-RLK, AGL79, and TCTP) were identified, which might probably be related to the formation of MMS phenotype. The present results not only would give a new insight into the genetic basis underlying the regulation of MMS but also would provide clues for plant architecture breeding in rapeseed.


INTRODUCTION
Plant architecture is a significant agronomic trait that indirectly influence the seed yield by affecting the formation of photosynthetic product and storage of grain-filling substance. Notably, the first "Green Revolution" was one of the great successes, which had improved plant architecture, especially in significant reduction plant height and an increase of the overall seed yield and harvest index (Khush, 2001). Seed yield per unit area is closely related to many plant architecture elements, including yield component, plant height, tiller number, intersection angle of tiller, and leaf shape. Especially in rice and wheat, the increase of the effective tiller number could improve their seed yield per plant and thus increase the overall yield. Some studies that are focused on the mechanism of plant architecture formation have been reported in rice, Arabidopsis thaliana (A. thaliana) and tomato (Leibfried et al., 2005;Jiang et al., 2013;Zhou et al., 2016). Therefore, improving the plant architecture of crops could effectively increase the crop yield.
Rapeseed (Brassica napus, AACC genome, 2n = 38) originated from a spontaneous hybridization between Brassica rapa (AA,2n = 20) and Brassica oleracea (CC, 2n = 18) (Nagaharu, 1935). Rapeseed is cultivated throughout the world for the production of vegetable oil, animal feed, and biodiesel. At present, the proportion of rapeseed oil in the total oil production of crop is as high as 57.2% in China (Xie et al., 2005;Yin et al., 2010); in the meanwhile, 60% of edible oil have to be imported from abroad (Wang and Yin, 2014). Seed yield of rapeseed is directly determined by yield-component traits, including seed number per silique, seed weight, and the number of effective siliques per plant (Qzer et al., 1999;Quarrie et al., 2006). Several indirectly related traits, including plant height, the number of the first effective branch, and the silique number of main inflorescence, were also important contributions. The genetic basis of the seed yield components and seed yield-related traits have been dissected by quantitative trait loci (QTL) mapping (Chen et al., 2007;Ding et al., 2012;Zhao et al., 2016) and genome-wide association studies (GWAS) (Liu et al., 2016;Lu et al., 2016).
In recent years, the new agronomic traits, multi-main stem trait (MMS, including double-main stems, three main stems, and over three main stems), have been discovered in rapeseed. MMS is an important seed yield-related trait, which was differentiated from shoot meristem. The plants with MMS have many advantages, such as multiple main stems (or tiller number), higher growth potential, fewer branches, and more seed number, which play an important role for increasing seed yield in rice and wheat (Lu et al., 2015;Shao et al., 2019;Hu et al., 2017). Favorable main stem number can also increase planting density, which contributes to a higher yield per unit area. In the process of differentiation, the morphological structure, physiological status, and endogenous hormonal content of the shoot apical meristem (SAM) undergo substantial change (Gordon et al., 2009). The SAM must maintain a balance between stem cell niches updating and peripheral organ initiation in order to fulfill this function. The genes that cause MMS might have a relationship with WUSCHEL (WUS) and CLAVATA (CLV) (Brand et al., 2002;Müller et al., 2006). In normal condition, the gene WUS is always expressed in a small group of cells in the organization center of the SAM, underneath the three outer cell layers (Weigel and Jürgens, 2002;Williams and Fletcher, 2005). Cells in center zone proliferate to form two kinds of cells (Haecker and Laux, 2001): one is progeny of stem cell holding multipotency, which maintain the structure of center zone, and the other is a daughter cell, who will develop into perimeter zone and then differentiate into various organizations. Only when the two progresses develop at the same pace, SAM will develop into normal structure; otherwise, the SAM will abnormally develop . During development process of SAM, WUS functions in maintaining the stem cells undifferentiated to keep natural structure and function of SAM (Schoof et al., 2000). Loss of function of WUS would lead to a large amount of stem cells developed into perimeter zone, and then giving rise to the leaf primordiums and secondary SAMs that will come into being . Another important gene in SAM is CLV that can promote stem cells to develop into perimeter zone. A WUS-CLV feedback loop in plant exists, which is of great importance to maintain normal structure of SAM. CLV3 encodes a kind of secreted protein to activate CLV1/CLV2 complex when the number of stem cells in SAM increase. When WUS is inhibited, the number of stem cells will decrease which in turn restrains the expression of CLV3 Dodsworth, 2009). Based on the previous reports, the maintenance homeostasis of SAM mainly includes WUS-CLV feedback loop (Fletche, 2018), phytohormone (Brand et al., 2000), redox regulation (Schippers et al., 2016), and network consisting of various genes.
QTL mapping has been proved to be an important approach for dissecting complex traits (Paran and Zamir, 2003). Many agricultural traits, for example, seed yield (Zhao et al., 2016), seed yield-related traits (Shi et al., 2009;Zhao et al., 2016), oil content , and fatty acid composition (Jiang et al., 2014;Bao et al., 2018), have been dissected by QTL mapping. In addition, some important genes that regulated agricultural traits of rapeseed were identified by fine mapping and map-based cloning, such as recessive genic male serility BNAMS1 (Yi et al., 2006), seed weight TSWA5a and TSWA5c (Fan et al., 2010), and seed number per silique BnaC9.SMG7b , etc. Gene-Fishing technology, based on PCR technology, could identify the differentially expressed genes quickly (Kim et al., 2004). For example, Choi et al. (2008) found 31 differentially expressed fragments in the grapevines inoculated with Rhizobium and salicylic acid by using Gene-Fishing (Choi et al., 2008). Junaedi et al. (2007) explored the effects of high allelopathic rice varieties on barnyard grass and discovered nine differentially expressed fragments (Junaedi et al., 2007). Park et al. (2006) screened the related genes that regulated seed epidermal needling of carrot, and 11 differentially expressed fragments were identified (Park et al., 2006). In addition, some other genes, including HIP1, CSH1, PPTSPO1, and SLFN-2, were also identified by Gene-Fishing technique Bradley et al., 2007;Frank et al., 2007;Sohn et al., 2007).
In the present study, the MMS characteristics were investigated in a doubled-haploid population named as KN DH population. We counted the number of multi-main stem in each DH line in Yangling and Wuhan for two successive years. Based on QTL mapping and Gene-Fishing technique, the purposes of the present study focused on two aspects: (1) identify QTLs for MMS and obtain candidate genes within QTL regions (2) based on differentially expressed genes and candidate gene within QTL regions and identify several common candidate genes that regulated the MMS.

Plant Materials and Field Experimental Design in Two Different Ecological Regions
The KN doubled-haploid (DH) population were derived from bi-parental segregating populations by microspore culture (Figure 1) . The KN population and its parents were sown in the middle of September in the experimental field of Hybrid Rape Research Center of Shaanxi Province, Yangling of Shaanxi province (YL) in northwest China (34°16'N, 108°5'E), and in the early October in the experimental field of Huazhong Agriculture University, Wuhan of Hubei province (WH) in central China (30°47'N, 114°35'E), for two consecutive years from the years 2014 to 2015, respectively. Yangling belonged to the winter ecological region, and Wuhan belonged to the semi-winter ecological region. The KN population and their two parents were planted in Yangling and Wuhan with two replications, respectively. Meanwhile, each particular combination of experimental year × location was defined as an independent experiment. Each independent experiment consisted of 348 DH lines. Each line was grown in a two-row plot with 40 cm between rows and 20 cm between individuals, and row length of 250 cm in all independent experiments. The field experiments followed a randomized complete block design and the normal agricultural practice.

Phenotypic Statistics and Analysis of MMS
The total number of each line corresponding to the number of plant with MMS in each line was investigated, and the ratio of the MMS was calculated in Yangling of Shaanxi province and Wuhan of Hubei province for two consecutive years. The ratio of the MMS was calculated as the ratio of the number of the MMS to the total number of each line. Not all of the DH lines had MMS in the KN population. In addition, the variance analysis and broad-sense heritability (h2) was calculated by using SAS 8 (SAS Institute Inc, 2000). The formula was h2 = б2g/(б2g + б2ge/n+ б2e/nr) × 100%, б2g is the variance among DH lines, б2ge is the interaction variance of the genotype with environment, б2e is the error variance, n is the number of environments, and r is the number of replications. Meanwhile, a DH lines with double-main stem (DMS) in KN DH population were selected for germination and morphological observation. The early morphology was observed in the five growth stages of seedling stage (15, 25, 30, 35, and 50 days). The seedlings for early morphological observation of the DMS materials and their natural plant type were grown in the culture chamber with 16 h light per day at 22°C ± 4°C.

QTL Mapping and Identification of Candidate Genes
The ratio of MMS in the KN population was collected as phenotype and used for QTL mapping. A high-density genetic linkage map of the KN population with 3,207 markers, including single-nucleotide polymorphism (SNP), SSR (simple sequence repeat), STS (sequence-tagged site), SRAP (sequence-related amplified polymorphism), and IFLP (intron fragment length polymorphism) markers, was constructed. The total length reached 3072.7 cM, and the average distance was only 0.96 cM between adjacent markers (Chao et al., 2017). Combining the ratio data of MMS with the high-density linkage map of KN population, QTL detection for MMS was performed using a composite interval mapping (CIM) with Windows QTL Cartographer 2.5 software (Wang et al., 2007). The scan walking speed was 1 cM; the window size was 10 cM with five background cofactors. According to the 1,000-permutation test corresponding to P = 0.05, the LOD value 2.0 was used as threshold for detection of significant QTLs. All detected QTLs were denoted as significant identified QTLs (Burns et al., 2003). The method for QTL nomenclature was as described by Wang et al. (2013). QTL integration was by metaanalysis using BioMercator 4.1 (Arcade et al., 2004). One QTL that was detected in at least two environments was taken to be a stable QTL, otherwise a specific QTL.
The alignment of the genetic map to the physical map and the identification of candidate genes were the same as those described by Chao et al. (2017). According to the collinearity of the highdensity genetic map and B. napus "Darmor-bzh" reference genome (Chalhoub et al., 2014), genome regions corresponding to the QTL confidence intervals (CIs) were identified by using closely linked SNP within QTL CIs. The candidate genes within QTL CIs were regarded as candidate genes (Shi et al., 2015). The orthologs of A. thaliana involving SAM were gathered from previous reported . The orthologous of candidate genes and their annotation were obtained by BLASTn based on A. thaliana database (http://www.arabidopsis.org/).

Additive × Additive Epistatic Interactions for MMS
To estimate the epistasis interaction of MMS, we used inclusive ICIM method with IciMapping software to analyze epistatic loci pairs Meng et al., 2015). To identify significant epistasis interaction loci, a walking speed was set to 1 cM, and significant LOD threshold was set at the same as default parameter of 5.0. The epistasis interaction loci were defined by abbreviation "Ep" of epistasis and trait name and chromosome number as well as location on chromosome (e.g., EpMMS.A8-115).

Acquisition of Differential Expression Genes Based on Gene-Fishing™ Technology
Differential expression genes (DEGs) were identified by using Gene-Fishing ™ DEG Premix Kits (Seegene, Seoul, South Korea), with an annealing control primer (ACP)-based PCR method (Kim et al., 2004). The total RNA of the DMSs was extracted and then reversed transcription for cDNA by using ReverTra Acea-(code no. FSK-100, TOYOBO). According to Gene-Fishing ™ PCR liquid system (20µl total volume system included cDNA 5µl, arbitrary ACP [5 µM] 2µl, dT-ACP2 [10 µM] 1µl, distilled water 2µl, 2×SeeAmp ™ ACP ™ Master Mix 10µl). The amplified PCR products were separated by 2% agarose gels electrophoresis with ethidium bromide, and differential expression fragments were discriminated. Differential expression fragments were isolated by using AxyPrep DNA Gel Extraction Kit (TIANGEN) and connected to a pMD 18-T Simple Vector (TaKaRa) and were sequenced. Sequencing data were confirmed with the GenBank database through the BlastX program of NCBI (http://www. ncbi.nlm.nih.gov/BLAST/) and/or the A. thaliana and B. napus database. Then, the sequence differences of genes were compared in natural plant and the DMS plant.

Quantitative Real-Time PCR (qPCR) Analysis of Six Potential Candidate Genes Regulating MMS
In order to investigate the potential candidate genes regulating MMS, several common candidate genes were selected for quantitative real-time PCR analysis. The shoot meristem of the natural plant and DMS plant in 5, 10, 15, and 20 days after germination was collected and used for relative expression analysis. The total RNA extraction and reversed transcription were the same as above method. Gene-specific primers were designed by using Primer 5.0, and the primers of the six target genes and reference gene actin were listed in Table S1. The RT-PCR reaction system was performed with three technical replicates by using TOYOBO SYBR ® R Green Realtime PCR Master Mix (code no. QPK-201) Kit. The amplification program was as follow: 95°C for 10 min, then 40 cycles with 95°C for 15s, 60°C for 15 s, and 72°C for 30 s, the last 60 to 95°C to do the melting curve.

Phenotype Investigation of KN DH Population and Performance of Multiple Main Stem
We investigated the plant with MMS for each line in the KN DH population for two consecutive years from the years 2014 to 2015 in Wuhan and Yangling. According to investigations in the years 2014 and 2015, 45 and 76 DH lines showed MMS phenotype in Yangling; however, only 25 and 29 DH lines showed MMS phenotype in Wuhan (Figure 2, Table 1). Only three lines (QT22, QT229, and QT236) of the 348 lines were found with MMS in all four environments. According to the number of plant with MMS and the total number of each line, the ratio of the MMS occurrence was calculated, the highest ratio of MMS was 100% in Yangling in the years 2014 and 2015, and in Wuhan in the year 2015, the lowest ratio was 20.1% in the year 2015 in Yangling (Table 1). In addition, 7.18-21.84% of plants in KN population had the phenotypes of MMS (from two to three main stems), with more leaves in comparison to their parents (Figure 1, Table 1). Furthermore, broad-sense heritability (h 2 ) was also calculated for MMS (Table S2). MMS had an h 2 value 62.26% in Yangling, and an h 2 value of 70.23% in Wuhan; these results indicated that the MMS had greater sensitivity and plasticity and was easily influenced by environmental condition.
Moreover, the plant with DMS and the natural normal plants from QT22 were used for germination and morphological development observation in 15, 25, 30, 35, and 50 days of different seedling stages (Figure 3). It was revealed that from germination to 25 days, no obvious differences were observed between the DMS plant and the normal plant, and DMS was first observed in 30 days. The stem base of the DMS plant hypertrophied gradually and had formed two main terminal buds. In 35 days, the DMS plant had formed two stems and could be obviously observed. In 50 days, the two main stems of DMS plants developed simultaneously and formed a plant type of DMS. These results indicated that the genes regulating the formation of the DMS function was in early developmental stages. The DMS plants had also grown more leaves and fewer branches compared to the natural plant (Figure 1, Figure 3E2).

QTL Mapping for MMS in the KN Population and Candidate Genes Identification
The MMS ratio in Yangling and Wuhan of the years 2014 and 2015 was used for QTL identification. Totally, 43 identified QTLs for MMS were detected in the KN DH population, which were located on 14 chromosomes, except for A02, A09, C02, C04, and C08 ( Figure 4, Table 2, Figure S2-1, Figure S2-2). It was revealed that q14WH7-2, q14WH15-1, and q15WH13-2 could explain the phenotypic variation more than 10%, which reached 11.38, 11.42, and 14.9%, respectively. Otherwise, QTL q14WH10-1 and q14WH10-2 were with the smallest phenotypic variation of 2.95%. In addition, the additive effects of QTLs were from −0.14 to 0.09, and the average CI was 4.77 cM ( Table 2). QTLs for MMS with overlapping CIs were integrated into consensus QTLs by using BioMercator 4.2 software. Forty-three identified QTLs were integrated into 41 consensus QTLs, of which 39 consensus QTLs that were only expressed in Yangling or Wuhan, were considered as environmental specific QTLs (Table 2, Figure 4, Figure 5). The remaining two consensus QTLs (cqMMS.A3-2 and cqMMS.C3-5) were repeatedly detected in two successive years in Yangling and Wuhan (Table 2, Figure 5), respectively, which were considered to be environmental stable QTLs. These results indicated that the majority of consensus QTLs for MMS were expressed in a specific environment, and MMS was greatly affected by environmental conditions. Regretfully, no major QTL for MMS was identified in all environments.

Additive × Additive Epistatic Interactions for MMS in Different Chromosomes
To estimate the epistasis interaction for MMS, epistatic loci pairs were detected by using IciMapping software. Totally, 26 statistically significant epistatic loci pairs were detected in three plant regions (Figure 6,  , had positive epistatic effect, indicating that the effects of parents for MMS were larger than their the recombinant effects of epistatic interaction pairs. The remaining epistatic interaction pairs had negative epistatic effects, indicating that the recombinant effects of these epistatic interaction pairs for MMS were larger than their parental effects. Furthermore, we noted that seven epistatic interaction loci were located in the above-mentioned five QTL intervals. EpMMS.A7-105 was located in QTL region of cqMMS.A7-3, indicating that cqMMS.A7-3 also interacted

Differentially Expressed Genes Between the Natural Plant and the DMS Plant by Gene-Fishing™ Technique
Total RNA of main stem apex of the natural plant and DMS in seedling of 20 days were extracted by using Invitrogen TRIzol Kit. The dT-ACP1 primers from Gene-Fishing technique were used as primers for reverse transcription. In order to identify differentially expressed genes (DEGs) in DMS and the natural plant of B. napus, the PCR amplification of cDNA using 36 random amplified primers was conducted. Totally, 46 differential expression bands were found, including 28 down-regulated and 18 up-regulated in the DMS compared to the natural plant (Figure 7, and a full scan of the entire original gels attached as Figures S1-1-Figure S1-6; Table S5). According to the sequence alignment results, 31 genes were obtained (Table 4). Based on functional annotation, seven underlying genes (DEG1, DEG2, DEG4, DEG5, DEG34, DEG36, and DEG46) were annotated to regulate the apical meristem, which might be related to the formation of DMS phenotypes (Table 4). In addition, some genes that are related to biotic stimulus, defense response, drought, and salt-stress responses, as well as cold-response functional genes were also identified ( Table 4). Subsequently, genes that regulate apical meristem were conducted for further analysis. DEG1 is homologous to the Arabidopsis gene AT4G29040 (RPT2A), and it had five copies distributed on A01 (BnaA01g07860D), A08 (BnaA08g13550D), C01 (BnaC01g09460D), C07 (BnaC07g41720D), and C08 (BnaC08g13300D). DEG2 had two copies found on A10 (BnaA10g29350D) and C09 (BnaC09g33500D) and are homologous to the Arabidopsis gene AT5G57950 (halted root, HLR). DEG4 (BnaA06g24200D) is homologous to the Arabidopsis gene AT4G01710 (crooked, CRK) and encodes ARP2/3 complex 16-kDa subunit. DEG5 had one copy on A01 (BnaA01g02910D), which is homologous to the Arabidopsis gene AT4G34220, and it was considered to be a leucine-rich repeat receptor-like protein kinase gene (LRR-RLK). DEG34 had two copies located on A06 (BnaA06g40690D) and C07 (BnaC07g47290D), and they were homologous to the Arabidopsis gene AT4G39270, which has the same function with DEG5 and was also considered to be leucinerich repeat receptor-like protein kinase family gene (LRR-RLK). DEG36 is homologous to the Arabidopsis gene AT3G30260 (AGAMOUS-like MADS-box protein 79, AGL79), which has five copies found on A06 (BnaA06g30700D), A09 (BnaA09g02740D), C02 (BnaC02g38130D), C07 (BnaC07g25980D), and C09 (BnaC09g02190D). DEG46 (BnaC03g39720D) is homologous to the Arabidopsis gene AT3G16640 (TCTP), which is mainly expressed in meristematic regions of the shoot and root.

Relative Expression Level of Candidate Genes Based on QTL Mapping and Gene-Fishing
Based on the common candidate genes that were obtained from QTL mapping and Gene-Fishing technique, six commonly identified genes were annotated to be related to the formation of MMS, including RPT2A, HLR, CRK, LRR-RLK, AGL79, and TCTP, and these genes were further used to perform the relative expression analysis by quantitative real-time PCR between the normal plant and the DMS plant in four stages of seedling (5, 10, 15, and 20 days of seedling age) (Figure 8). It was revealed that the expression level of RPT2A was extremely and significantly higher in DMS plant than that of in natural plant in 15 days of seedling age, and in 10 and 20 days of seedling age, the expression level of RPT2A in DMS plant was extremely and significantly lower than that of in natural plant. For HLR, the expression level in DMS plant was much higher than that in natural plant in all four states of seedling age. Especially in 5, 10, and 15 days, the expression levels in DMS plant were significant or extremely and significantly higher than that in natural plant. In 10 to 20 days of FIGURE 5 | Environmental stable QTL and environmental specific QTL on A03, A07, and A08 chromosomes. In the left, orange QTL cqMMSa3-2 represents environmental stable QTL in Wuhan and Yangling. Middle green and red QTLs represent specific QTL in Wuhan. Right yellow QTLs represent specific QTL in Yangling. September 2019 | Volume 10 | Article 1152 Frontiers in Plant Science | www.frontiersin.org   Continued) seedling age, expression level were higher in natural plant than that in DMS plant, which were also consistent with Gene-Fishing results ( Table 5). In addition, whether in natural plant or the DMS plant, expression of HLR also showed downward trend in all four states. For CRK, in 5 days of seedling age, expression level in DMS plant was extremely and significantly higher than that in natural plant; however, in 10 days of seedling age, the expression of CRK was severely inhibited in DMS plant. From 15 to 20 days, expressions of CRK were slowly rising again. The decline in gene expression of CRK might lead to a decrease in the activity of some regulated genes, which disrupted the balance of SAM.
For LRR-RLK, the expression level in DMS plant was extremely and significantly lower than in natural plant of 10 and 20 days of seedling age; however, its expression level in DMS plant was extremely and significantly higher than in natural plant in 5 and 15 days of seedling age. In addition, the expression level of LRR-RLK in natural plant was the highest at fourth stage, which was compatible with result of Gene-Fishing (Table 5). For AGL79, the expression level in DMS plant was extremely and significantly lower than in natural plant in 5 days of seedling age. With the development of seedlings, the expression level in DMS plant was slightly higher (not significant) than in natural plant of 10 to 20 days of seedling age. For TCTP, the expression level in DMS was extremely and significantly lower than in natural plant of 5 to 15 days of seedling age. In 20 days of seedling age, the expression level of TCTP showed upward trend in DMS plant. Relative expression level of candidate genes showed that these six genes had different expression profile between the natural plant and the DMS plant in different seedling age stages.

DISCUSSION
Rapeseed is an important oilseed crops worldwide (Cai et al., 2014). Although oil content and seed yield of rapeseed varieties have been effectively improved by rapeseed breeders (Zhao et al., 2011), the optimal plant architectural cultivated varieties are urgently needed for mechanized harvesting and high-density planting in the current agricultural management (Fu, 2008). Rapeseed plant architecture is characterized by plant height, branch number and angles, and inflorescence morphology (Li et al., 2018a), which indirectly influence rapeseed yield by affecting the major yield-component trait (Qiu et al., 2006;Zhou et al., 2014). Crop breeding has primarily focused on plant architecture, including plant height, branch or tiller number and angle, leaf shape, and size . MMS differentiated from shoot meristem, which is also one of the important factors for improving plant architecture and is conducive to increasing plant density and enhancing seed yield in rapeseed. In KN DH population, several MMS plants were found, and their main characteristics were composed of double-or three main stems from the base of a plant and had less lateral branches (Figure 1). Because of nutrient diversion, the MMS plant had slight lower plant height and stronger stem; thus, the MMS materials had stronger lodging resistance. In addition, under a small sowing amount of seeds, the multiple main stem materials not only can effectively increase plant density and the total number of siliques per unit, which finally might increase seed yield, but also, it is suitable for mechanized harvesting and to improve production efficiency. In total, the acquisition of the MMS material provided completely new ideas to improve the plant architecture and breed new varieties in crops.
According to the statistics of the number for MMS materials in different plant environments in several years, MMS was a complex quantitative trait and was significantly influenced by the environment. QTL mapping is considered to be an effective approach for dissecting quantitative traits (Paran and Zamir, 2003). In the present study, 43 QTLs for MMS that were located on 14 chromosomes were identified, explaining 2.95-14.9% of the phenotypic variation ( Table 2). Two environment-stable QTLs, including qMMS.A3-2 and qMMS.C3-5, were expressed in both winter and semi-winter environments. More importantly, according to epistatic interaction analysis, QTL cqMMS.C3-5 had three epistatic interactions with other loci, which indicated that this stable QTL was an important loci for MMS. Several environmentstable QTLs for other agricultural traits in rapeseed have been reported, such as stable QTLs for flowering time (Li et al., 2018b;Shi et al., 2009) and thousand seed weight (Zhao et al., 2016). These environment-stable QTLs are suitable for developing new varieties of rapeseed with broad adaptability (Li et al., 2018b). Some genes within QTLs related to cell division and differentiation of shoot meristem were identified mainly included STM, WUS-related homeobox genes, cytokinin response factors, and auxin-related genes. Previous studies showed that cytokinin (CK) produced by root can promote cell division as well as maintenance of meristem, and the upward transport of CK in the root promotes axillary bud germination and the formation of lateral branches (Chatfield et al., 2000;Leyser, 2003;Ongaro and Leyser, 2008). Gene-fishing technology has been widely used to identify DEGs in crops (Sanghoon et al., 2009;Liao et al., 2015). Some important differently expressed genes have been identified, such as salt-stress-induced up-regulated genes in barley leaves (Sanghoon et al., 2009) and pistillody genes in wheat (Liao et al., 2015). In the present study, 31 differentially expressed genes between DMS plants and the natural plants of B. napus were identified, in which six DEGs involved in the formation and regulation of DMS phenotype (Table 4), including RPT2A, HLR, CRK, LRR-RLK, AGL79 (AGAMOUS-LIKE 79), and TCTP (translationally controlled tumor protein). Among them, RPT2A is a 26S proteasome subunit and was located in QTL regions of cqMMS.A1-3, cqMMS.  and was reported to be involved in maintenance of the postembryonic meristems (Ueda et al., 2004). 26S proteasome is an important factor in signal transduction-mediated pathway; missing or mutation of 26S proteasome would lead to a decrease in strigolactones (SL) content, thereby promoting the excessive growth of axillary buds and other organs, such as branches and tillers (Wang and Bouwmeester, 2018). Therefore, changes of 26S proteasome may affect the growth of SAM and result in MMS phenotype. Abnormal expression of RPT2A in DMS plant might change the state of meristem and SL content in 5 to 20 days of seedling age, thereby forming MMS phenotype. HLR is located in QTL regions of cqMMS.A10-1 and cqMMS.C9-1, which encodes 26S proteasome regulatory subunit RPT2A, which is essential to maintain the size of SAM and might affect the growth of SAM by affecting the expression of WUS (Ueda et al., 2004). However, the expression level of HLR in DMS plant was much higher than that in natural plant in four states of seedling age, and overexpression of HLR might disturb the normal state of meristem and initiate bud differentiation. Meanwhile, HLR is also required to regulate the cell division pattern in SAM. In plants, the actin cytoskeleton is essential for various processes such as stomatal closure, cell proliferation, and cell morphogenesis (Akin and Zipursky, 2014). CRK encodes ARP2/3 complex 16-kDa subunit. In moss Physcomitrella, when the constituent part arpc4 of the CRK is lost, shoot growth of the apical cell slows down obviously. In the current study, CRK was located in QTL region of cqMMS.A6-1, and the decline in gene expression of CRK might lead to a decrease in the activity of some regulate genes in seedling. Studies in A. thaliana, when the CRK is absent, the actin of trichomes, hypocotyls cells, and epidermal pavement cells gather into bundles, downstream cell activity was interfered (Deeks and Hussey, 2003;Quatrano et al., 2006), which disrupted the balance of SAM. LRR-RLK represented a large family of protein kinases; it was located in QTL region of cqMMS.A1-2 (Gou et al., 2010). Expression of LRR-RLK had upward trend in DMS plant, which indicated that LRR-RLK might participate in the occurrence of axillary buds (Yaginuma et al., 2011). In addition, LRR-RLK contains ERECTA (ER) and CLV1. ER-family receptor kinases can regulate stem cell homeostasis to regulate organ shape and inflorescence architecture via cushioning its CK response in the SAM (Shchennikova et al., 2017). CLV1 interacts with the WUS gene to maintain the balance of differentiation of shoot and floral meristem cells in SAM (Gou et al., 2010). When CK stimulates the expression of WUS, the expression of CLV3 should have been inhibited to result in a SAM expansion; however, the FIGURE 7 | Differential expressed genes screen using ACPs primers based on Gene-Fishing. The numbers and arrows on image represent different differential expression genes. For better presentation, only cropped images of blots are shown here; a full scan of the entire original gels is submitted in the Supplementary Material (Figures S1-1-S1-6). September 2019 | Volume 10 | Article 1152 Frontiers in Plant Science | www.frontiersin.org  morphology of SAM showed no difference with the normal plant. These phenomena could owe to the cushioning mechanism of ER family, which prevents the expression of CLV3 from increasing remarkably, thereby maintaining the normal development of stem meristem (Torii et al., 1996;Uchida et al., 2013). AGL79 is a transcription factor; relatively lower AGL79 overexpression would promote more lateral branches and rosette leaves (Gao et al., 2018). In addition, miR156/SPL10 modulated lateral root development and branching and leaf morphology in A. thaliana by silencing AGL79 (Gao et al., 2018). TCTP is a central regulator of cell proliferation and differentiation in animals and in plants, which is an important mitotic regulator (Amson et al., 2013;Roberto et al., 2015). Furthermore, some genes also have similar functions in other crops. HLR gene is a homolog of the 26S proteasome subunit, which was strongly expressed both in the SAM and RAM of rice seedlings (Yanagawa et al., 2002). LRR-RLKs play an important role in regulating the SAM and microspores (Becraft, 2002), brassinosteroid perception, and floral abscission (Carles and Fletcher, 2003;Li, 2003). Rameneni et al. performed the wholegenome sequencing (WGS) of B. rapa and identified the LRR-RLK of B. rapa, which is involved in the plant morphological characters and plant stress resistance (Rameneni et al., 2015). ERECTA is one of LRR-RLK family, which coordinates Arabidopsis organ growth and flower development by promoting cell proliferation (Torii et al., 1996;Shpak et al., 2004). VEGETATIVE1 (VEG1) is an AGL79-like MADs-box gene; the secondary inflorescences in nutation were replaced by vegetative branches in pea (Berbel   et al., 2012). In addition, AGL79 is involved in regulating leaf shape, shoot branching, and flowering time in A. thaliana (Gao et al., 2018). In addition, when maize is under stress response, TCTP expression showed a significant upward trend (Chen et al., 2014). Taken together, a series of balance will be broken when these genes mutated and may lead to disturbance in the growth of SAM and form abnormal organ morphology, such as clumping leaves and multiple axillary buds (Gou et al., 2010;Uchida et al., 2013).

CONCLUSION
The increase of yield and oil content is the ultimate goal in rapeseed breeding. Due to the lack of mechanized cultivated variety and high productive cost as well as low income, farmers prefer to idle the farmland rather than plant rapeseed, which seriously hampers the development of rapeseed in China. Importantly, the acquisition of MMS materials provides the possibility for breeding mechanized varieties, high-density planting varieties, and the yield improvement in rapeseed. In the present study, six common candidate genes within QTL regions and differentially expressed genes, including RPT2A, HLR, CRK, LRR-RLK, AGL79, and TCTP were identified by using QTL mapping and Gene-Fishing technique. Meanwhile, according to their functional annotation, these candidate genes might be related to the formation of MMS phenotype, in which genes might disrupt cell division and the differentiation of SAM and induce the formation of multiple buds, thereby promoting the formation and development of the MMS plant. MMS is an important agriculture trait, which is differentiated from shoot apical meristem. Relative expression level of candidate genes showed that these genes had different expression level between the natural plants and the DMS plants in different seedling age stages. Overall, the study would not only give a new insight into the genetic basis underlying the control of the MMS in rapeseed but also provide clues for plant architecture breeding in rapeseed.

DATA AVAILABILITY
All datasets generated for this study are included in the manuscript/Supplementary Files.

ETHICS STATEMENT
The authors declare that the experiments comply with the current laws of the country in which they were performed.

AUTHOR CONTRIBUTIONS
WZ wrote the first draft of the manuscript. WZ and HC conducted all the field experiments, performed most of the experiments and data analysis for the overall study; HC, LnZ, NT, YZ, BL, KZ, ZG and HW helped and participated in the field trials and data collection for multiple-main stem. KC, LbZ and DH performed the phenotypic data processing and QTL detection. ML designed and conceived the overall study and revised the manuscript. All authors read and approved the final manuscript.