Skip to main content

Identification of genomic regions associated with multi-silique trait in Brassica napus

Abstract

Background

Although rapeseed (Brassica napus L.) mutant forming multiple siliques was morphologically described and considered to increase the silique number per plant, an important agronomic trait in this crop, the molecular mechanism underlying this beneficial trait remains unclear. Here, we combined bulked-segregant analysis (BSA) and whole genome re-sequencing (WGR) to map the genomic regions responsible for the multi-silique trait using two pools of DNA from the near-isogenic lines (NILs) zws-ms (multi-silique) and zws-217 (single-silique). We used the Euclidean Distance (ED) to identify genomic regions associated with this trait based on both SNPs and InDels. We also conducted transcriptome sequencing to identify differentially expressed genes (DEGs) between zws-ms and zws-217.

Results

Genetic analysis using the ED algorithm identified three SNP- and two InDel-associated regions for the multi-silique trait. Two highly overlapped parts of the SNP- and InDel-associated regions were identified as important intersecting regions, which are located on chromosomes A09 and C08, respectively, including 2044 genes in 10.20-MB length totally. Transcriptome sequencing revealed 129 DEGs between zws-ms and zws-217 in buds, including 39 DEGs located in the two abovementioned associated regions. We identified candidate genes involved in multi-silique formation in rapeseed based on the results of functional annotation.

Conclusions

This study identified the genomic regions and candidate genes related to the multi-silique trait in rapeseed.

Background

In flowering plants, the pistil develops into a fruit, a process crucial for both yield and seed quality. Many recent studies have focused on abnormal pistil development, such as the formation of pistillodes (sterile pistils), pistil number variation, and mutations affecting pistil length [1,2,3,4,5,6]. Of these, pistil number variation, or the multi-pistil phenomenon, has also been reported in other plants, such as alfalfa (Medicago sativa) [1, 7], wheat (Triticum aestivum) [8, 9], sweet cherry (Prunus avium) [10], and rye (Secale cereale) [11].

Brassica crops are important sources of edible oil, vegetables, industrial erucic acid, animal feed (seed meal), and ornamental plants. Like most Brassicaceae family members, Brassica species have monoclinous flowers, with one pistil and six stamens per normal flower.

Morphological variation in siliques (pods) has been reported in several Brassica species. In the 1940s, Pathak [12] reported an abnormal plant of Brassica campestris L. var. sarson prain, which had one to two extra pods telescoped into the outermost pods. The inner pods was basically normal but were contorted and smaller than those of the wild type. Multi-locular Brassica juncea [13], trilocular B. juncea [14, 15] and Brassica rapa [16] were obtained from interspecific hybridization. Unlike previously reported fruit variations, zws-ms flower has three completely independent and apocarpous pistils and nine to ten stamens sharing the same receptacle in each floral organ, resulting in the formation of three independent siliques sharing the same carpopodium. Thus, zws-ms has an increased number of pod on each carpopodium rather than an increased number of loculi per pod. Similar phenotypes have also been observed previously: Guan and Li [17] obtained a line with double siliques through treatment of B. napus with 12C heavy ion beam irradiation. Hu et al. [18] identified a multi-pistil plant in a B. campestris L. cytoplasmic male sterile line, different from our zws-ms plants which were obtained by crossing of B. napus and B. rapa [19]. Though the agronomic and physiological traits, and genetic mechanism of zws-ms plants have been analyzed, the molecular mechanism underlying multi-pistil formation, which increased the silique number per plant in rapeseed, remains elusive.

Bulked segregant analysis (BSA) [20] is a simple but effective strategy for identifying molecular markers linked to target genes by genotyping a pair of bulked DNA samples from two populations with distinct phenotypes. BSA was usually combined with traditional markers, such as simple sequence repeat (SSR) [21], sequence related amplified polymorphism (SRAP) [22] or random amplified polymorphic DNA (RAPD) [23] markers to map traits. Recent years, combination of BSA and next-generation sequencing (NGS) have been successfully used in the identification of key genes/loci related to regular agronomic traits in B. napus. For example, Wang et al. combined BSA with NGS to map quantitative trait loci (QTL) for the branch angle and revealed BnaYUCCA6 was candidate gene controlling branch angle [24]. Geng et al. combined specific length amplified fragment sequencing (SLAF-seq) with BSA and identified four candidate genes related to seed weight [25]. Yao et al. mapped the orange petal color gene (Bnpc1) in spring B. napus to a 151-kb region via BSA with whole genome re-sequencing (WGR) [26]. Hence, we believed that combining BSA and WGR would be a convenient way to map the genomic regions responsible for the multi-silique trait.

In this study, we performed association analysis based on SNP and InDel data using a Euclidean Distance (ED) algorithm to identify the gene underlying the multi-silique trait in rapeseed. Though the SNP-index method is also a classical method for association region identification, and has been employed in many studies [25, 26]. However, it is suitable in the case when four DNA samples are sequenced: two parental plants and two extreme pools from segregated population. In this study, only two pools from NILs multi-silique pool and single-silique pool were sequenced, therefore, only ED algorithm was useful in this situation. We also performed transcriptome sequencing to identify differentially expressed genes (DEGs). Based on functional annotation of the genes with SNPs or InDels and DEGs in the associated regions, we identified candidate genes that might be closely related to the multi-silique trait.

Results

Morphological analysis of the multi-silique trait in zws-ms

Under normal growth conditions, the NILs (Fig. 1) zws-ms (multi-silique) and zws-217 (single silique) showed no phenotypic differences (e.g., plant height, leaf shape, and number of branches) at the vegetative stage. However, zws-ms had shorter and rounder buds than zws-217 (Fig. 2a). At the full-bloom stage, zws-ms showed distinct floral structures with zws-217. In contrast to normal Brassicaceae flowers, which typically contain a single pistil and six stamens (four tetradynamous stamens and two shorter stamens), zws-ms flowers had three completely independent apocarpous pistils sharing the same receptacle in each floral organ (Fig. 2b). The zws-ms flower had nine to ten stamens, six of which were more prominent than the others (Fig. 2b). In zws-ms flowers, all pistils and stamens shared the same receptacle, but the petals and calyxes appeared normal. Roughly 32~53% of flowers in a zws-ms plant had this mutant multi-pistil and stamen configuration, whereas the rest flowers in the same plant were normal. The abnormal zws-ms flowers developed into three independent siliques sharing the same carpopodium as the flowers (Fig. 2c). The triple siliques were normal in length, but occasionally one became degenerated and curved.

Fig. 1
figure1

Scheme for constructing the multi-silique population (zws-ms) and the corresponding single-silique NIL population (zws-217)

Fig. 2
figure2

Morphological differences between zws-217and zws-ms. a: at the budding stage, zws-ms has inflated, round buds; b: at the full-bloom stage, a normal flower has one pistil and six stamens, while zws-ms has three pistils and nine stamens per floral organ; c: a normal silique in zws-217 and the three-silique trait in zws-ms. These images were taken originally in 2017

Association analysis

High-throughput whole genome re-sequencing on the Illumina HiSeq platform produced a total of 89.56 Gb of raw data. The Q30 ratio averaged 90.52%. After alignment to the reference genome, the average coverage was 91.78% and the average sequencing depth was 40×. When the clean reads were mapped to the B. napus reference genome, each of them was aligned in both forward and reverse directions; thus, the number of total reads was twice of the clean reads: zws-ms and zws-217 contained 328,507,566 and 263,588,948 reads, with 89.89 and 89.10% properly mapped to the reference genome, respectively.

After filtration, 1,135,690 SNPs and 289,832 small InDels of high quality (HQ) were used for association analysis. We used the ED algorithm due to only two pools from NILs zws-ms and zws-217 were sequenced. According to the ED algorithm [27], the association threshold was 0.26. Based on this information, three SNP-associated regions were identified, including one on chromosome A09 (Fig. 3a) and two on chromosome C08 (Fig. 3b). One smaller associated region on the chromosome C08 that close to another larger one (Fig. 3b) was also detected. These regions had a size of 10.45 MB and contained 2095 genes. Similarly, two InDel-associated regions were identified on chromosomes A09 (Fig. 3c) and C08 (Fig. 3d) totaling 10.20 MB and containing 2044 genes, with an association threshold value of 0.30. Combining the SNP-associated regions and InDel-associated regions, two merged intersecting regions (Table 1, Additional file 1: Figure S1) were obtained on chromosomes A09 and C08, with 2044 genes in this 10.20 MB region, including 126 genes with non-synonymous mutations and 23 with frame-shift mutations.

Fig. 3
figure3

Associated regions calculated using the ED algorithm. Black solid lines represent the fitted values; red dotted lines represent the associated threshold. a: SNP-associated region on ChrA09; b: two SNP-associated regions on ChrC08; the smaller region is indicated by a red arrow; c: InDel-associated region on ChrA09; d: InDel-associated region on ChrC08

Table 1 Information about the two intersecting regions

Transcriptome sequencing, data mapping, and annotation

We performed sequencing saturation and cluster analysis of the samples to ensure the validity of the data. In total, 61 Gb of raw data were generated, with an average Q30 value of 90.54%. For each sample, approximately 36.7 M clean reads were generated, with an average GC content of 47.23% (Additional file 2: Table S1). The average proportion of total reads mapped to the reference genome of samples T01 to T06 was 73.729% (Additional file 3: Table S2), ranging from 77.81 to 76.03%.

Annotation of genes containing SNPs and InDels

To reveal the biological functions of the potential candidate genes, we performed GO and KEGG pathway enrichment analyses. Within the SNP- and InDel-associated regions, 2041 of 2044 genes could be annotated (Additional file 4: Table S3). GO annotations included 2217 terms involved in biological processes (BP, Fig. 4, Additional file 5: Table S4), 326 terms in the cellular component (CC, Fig. 4, Additional file 6: Table S5), and 799 terms in the molecular function (MF, Fig. 4, Additional file 7: Table S6). The BP terms with the highest levels of enrichment included response to nitrate (GO: 0010167), nitrate transport (GO: 0015706), negative regulation of flower development GO: 0009910), transition metal ion transport (GO: 0000041), inositol biosynthetic process (GO: 0006021), and organ development (GO: 0048513) (Additional file 5: Table S4). The CC category includes cell wall (GO: 0005618), vacuole (GO: 0005773), cytosolic small ribosomal subunit (GO: 0022627), and small ribosomal subunit (GO: 0015935) (Additional file 6: Table S5). In the MF category, the most enriched terms were sinapoyl spermidine: sinapoyl CoA N-acyltransferase activity (GO: 0080089), 2-isopropylmalate synthase activity (GO: 0003852), methyl salicylate esterase activity (GO: 0080031), nucleocytoplasmic transporter activity (GO: 0005487), and calmodulin binding (GO: 0005516) (Additional file 7: Table S6).

Fig. 4
figure4

Gene ontology (GO) terms of the genes from the intersecting regions. A total of 2041 genes were divided into three categories: biological processes, cellular components, and molecular functions

Moreover, 126 genes with non-synonymous mutations were identified, including 104 genes with GO annotations (Additional file 8: Table S7). Additionally, 23 genes with frame-shift mutations were identified, including 19 with GO annotations (Additional file 9: Table S8). Analysis of these annotations produced some interesting results: BnaA09g43100D (Bna.A09SFH3) and BnaC08g38270D (Bna.C08HAC12) were annotated as “flower development (GO:0009908)”; BnaC08g36330D (Bna.C08GPRI1) is related to “negative regulation of flower development (GO:0009910)”; BnaA09g44210D (Bna.A09BES1) and BnaA09g48960D are associated with “ovule development (GO:0048481)” and “carpel development (GO:0048440)”, respectively; BnaA09g42920D (Bna.A09EC1.3) was annotated as “regulation of double fertilization forming a zygote and endosperm (GO:0080155)”; and BnaA09g42700D was annotated as “specification of floral organ number (GO:0048833)”.

KEGG pathway enrichment analysis revealed the four most significantly enriched pathways involving the DEGs: stilbenoid, diarylheptanoid and gingerol biosynthesis (ko00945), limonene and pinene degradation (ko00903), glycosyl phosphatidyl inositol (GPI)-anchor biosynthesis (ko00563), and isoquinoline alkaloid biosynthesis (ko00950) (Additional file 10: Table S9). These pathways were further classified into five major groups: cellular processes, environmental information processing, genetic information processing, metabolism and organismal systems (Fig. 5). Of these, the sub-group plant hormone signal transduction had the highest number of annotated genes.

Fig. 5
figure5

Classified KEGG pathways of the genes from the intersecting regions. The pathways were further classified into five major groups: cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems

Analysis of differentially expressed genes and potential candidate gene identification

To identify DEGs related to the multi-silique trait, we examined the expression levels of transcripts from zws-ms and compared them with those of normal plants (zms-217). A total of 129 DEGs were identified in buds, with 52.94% (67 genes) upregulated and 48.06% (62 genes) downregulated (Additional file 11: Table S10), 39 of which (Table 2) were located in the two associated regions mentioned above. Thirteen genes were located on chromosome A09, while 26 were located on chromosome C08. These 39 genes were clustered into two groups based on expression level (Fig. 6, Table 2): the first group contained 14 genes that were upregulated in multi-silique zws-ms plants, while the second group contained 25 genes that were downregulated in these plants.

Table 2 GO annotations of the 39 DEGs located in the two associated regions
Fig. 6
figure6

Expression patterns of 39 DEGs located in the two associated regions. The expression levels are given in log2(FPKM+ 1). FPKM: Fragments per Kilobase Million. T01, T02, and T03: three random samples from zws-ms at the budding stage; T04, T05, and T06: three random samples from zws-217 at the budding stage

Of the 39 DEGs, three (BnaA09g45320D (Bna.A09CPN10), BnaC08g39130D (Bna.C08CPN10), and BnaC08g41780D (Bna.C08SRS)) were assigned to GO categories involved in ovule development (GO: 0048481), BnaC08g42450D (Bna.C08RLK7) in stamen development (GO: 0048443) and BnaC08g41720D (Bna.C08APA1) in organ development (GO: 0048513).

Validation of transcriptome sequencing data by qPCR

We validated the transcriptome sequencing data by comparing the relative transcript levels of these DEGs in zws-ms and zws-217 by qPCR. qPCR analysis of plants in buds showed that all genes except BnaC08g39120D had similar trends in expression as those observed by transcriptome sequencing, especially BnaA09g45320D and BnaC08g40410D (Figs. 6 and 7). These results further confirm the reliability of the transcriptome sequencing data described above.

Fig. 7
figure7

Relative expression levels (folds) of 10 selected genes determined by qPCR. Note: a: BnaA09g45320D; b: BnaC08g40410D; c: BnaA09g45890D; d: BnaC08g41720D; e: BnaC08g42080D; f: BnaC08g40740D; g: BnaA09g45310D; h: BnaA09g47900D; i: BnaC08g39120D; j: BnaC08g41780D

Orthologs of these 12 candidate genes in Arabidopsis

These 12 candidate genes mentioned above, including seven with non-synonymous mutations or frame-shift mutations, as well as the five DEGs, were then aligned with Arabidopsis and their orthologs (Table 3) in this well-studied model plant were identified based on The Arabidopsis Information Resource (TAIR, https://www.arabidopsis.org): (1) AT2G22370 (ortholog of BnaA09g42700D), identified as MED18, encodes a subunit of the mediator complex that affects flowering time and floral organ formation through FLOWERING LOCUS C (FLC) and AGAMOUS (AG); (2) AT2G21750, identified as EGG CELL 1.3 (EC1.3), encodes a small cysteine-rich protein secreted by the egg cell during gamete interactions; (3) AT2G21540 is an SEC14-LIKE 3 (ATSFH3) protein; (4) AT1G19350 (BZR2 or BES1) encodes brassinosteroid signalling protein that accumulates in the nucleus as dephosphorylated form in response to brassinosteroid. (5) BnaA09g45320D and BnaC08g39130D share the same ortholog, AT1G14980, known as CPN10, encoding mitochondrial-localized chaperonin 10; (6) AT1G08260, known as TIL1 gene, encodes the catalytic subunit of DNA polymerase ε (abo4–1); (7) AT2G20570 (GPRI1 or GLK1) encodes one of a pair of partially redundant nuclear transcription factors that regulate chloroplast development in a cell-autonomous manner; (8) AT1G11910 (APA1) encodes an aspartic proteinase that forms a heterodimer and is stable over a broad pH range; (9) AT1G11870 is identified as OVA7 (SRS), of which expression products Seryl-tRNA synthetase targeted to chloroplasts and mitochondria; (10) AT1G09970 (RLK7) belongs to a leucine-rich repeat class of receptor-likekinase (LRR-RLKs); (11) AT1G16710 (HAC12) encodes an enzyme with histone acetyltransferase activity that can use both H3 and H4 histones as substrates.

Table 3 The 14 candidate genes and their orthologs in Arabidopsis

Discussion

In the current study, we constructed NILs to establish the two extreme DNA pools (multi-silique pool and single-silique pool), instead of the frequently used F2 [28,29,30], doubled haploid (DH) [25, 31], or recombinant inbred line (RIL) [32, 33] populations. BSA [20] combined with NGS [34], which has been successfully used in various crops [29, 30, 32, 35, 36]. Here, the WGR produced abundant, high-quality data. After filtration, 1,135,690 HQ SNPs and 289,832 InDels were identified and used to conduct association analysis. These markers were clearly superior to traditional markers in terms of number, efficiency, and quality: as we expected, these markers were much denser than SSR or other markers and the combination of BSA and WGR saved the labors and shortened the whole process. As we were unable to use the SNP-index method in this study, we used the ED algorithm instead, followed by the identification of three SNP-associated regions, including one on chromosome A09 and two on chromosome C08. The two regions on chromosome C08 were very close to each other (only 0.14 Mb apart), and the smaller region was only 0.07 Mb long, suggesting that it is likely a minor site. Moreover, we identified two InDel-associated regions. Taken together, the SNP-associated region data were highly consistent in terms of position with those of the InDel-associated regions, confirming the reliability of both data sets. Therefore, the two intersecting regions are considered to be highly associated with the multi-silique trait, representing the only difference between the zws-ms and zws-217 pools. Preliminary analysis indicated that these two regions contain 2044 genes, including 126 with non-synonymous mutations and 23 with frame-shift mutations.

Since these genes were identified by WGR, they had small differences in sequence, which consequently led to changes in the amino acid sequences, structures, and final functions of the proteins. Since zws-ms and zws-217 are NILs that only differ from each other in the multi−/single-silique trait, we reasoned that these genes are associated with this trait and were therefore selected as candidate genes. However, some genes might have differed between the two pools in terms expression levels rather than nucleotide sequence. Therefore, we performed transcriptome sequencing from buds in plants to uncover the DEGs between lines, which we then confirmed by qPCR. We identified 129 DEGs between zws-ms and zws-217, which are distributed on all chromosomes except chromosome A06. Among the 129 DEGs, 39 are located in the two important trait-associated regions identified in this study and are therefore considered highly related to the multi-silique trait. These DEGs were classified into two groups based on expression level, including one containing upregulated genes and the other containing downregulated genes.

A comprehensive understanding of the functions of the 2044 genes in the two intersecting regions is crucial for candidate gene selection. Unlike standard agronomic traits (such as plant height and flowering time), which are relatively well-studied and for which some candidate genes are available for reference, the multi-silique trait analyzed in this study has not been studied on genomic or molecular level in detail to date. Thus, GO annotation and KEGG pathway analysis provided important clues for candidate gene identification. Among the 2044 genes within the two associated regions, 2041 were successfully annotated. KEGG pathway analysis provided a general view of the genes in these regions. Most enriched KEGG pathways were involved in the biosynthesis or degradation of some secondary metabolites.

As we were unable to relate this information to the multi-silique trait using established knowledge, we conducted a GO analysis. Eight of these DEGs were annotated in the category “negative regulation of flower development (GO: 0009910)” and 221 genes were involved in “organ development (GO: 0048513)”. This information is too general to reveal an association with flower/silique development. In addition, most genes within the two regions were the same between the two DNA pools. Thus, more attention should be paid to genes showing polymorphism between zws-ms and zws-217 plants, specifically genes harboring non-synonymous mutations or frame-shift mutations. Furthermore, 104 of the 126 genes with non-synonymous mutations and 19 of the 23 genes with frame-shift mutations were represented in the GO database. We examined each of their annotations and identified seven interesting genes, which were annotated to ovule, carpel or stamen development. The NILs zws-ms and zws-217 differ only in terms of flower and silique morphology and structure. Therefore, genes that differ in zws-ms versus zws-217 at the genomic level are thought to encode diverse proteins, suggesting they are excellent candidate genes for the multi-silique trait.

At the transcriptome level, we focused on DEGs, especially those located in the associated regions mentioned above. Thirty-nine out of 129 DEGs were found in the two intersecting regions, among which five genes were annotated to the ovule or stamen development. This annotation information strongly suggests that these DEGs are involved in the multi-silique trait, as zws-ms flowers have obviously abnormal morphology, consequently leading to the formation of abnormal fruits.

Then, we aligned these selected genes to the Arabidopsis, screening for homologous genes sharing the highest sequence identity with them, but did not identify more potential candidate genes other than these 12 genes already mentioned (BnaA09g43100D, BnaC08g38270D, BnaC08g36330D, BnaA09g44210D, BnaA09g48960D, BnaA09g42920D and BnaA09g42700D with non-synonymous mutations or frame-shift mutations, as well as BnaA09g45320D, BnaC08g39130D, BnaC08g41780D, BnaC08g42450D and BnaC08g41720D with different expression levels). Based on analysis of these orthologs in Arabidopsis and their roles in plants, the annotated information of these candidate genes are confirmed: (1) BnaA09g42700D is found much shorter than its Arabidopsis ortholog AT2G22370 (MED18) in length, but in their overlapping region, they share 100% identity of nucleotide sequence. Mutation of MED18 in Arabidopsis significantly delayed its flowering time [37]. FLC encodes an MADS-box protein and is a repressor for flowering [38]; AG is a floral meristem identity gene, which specifies carpel and stamen identity in the flower [39, 40]. Thus, BnaA09g42700D could probably act on flowering indirectly. (2) Disruption of AT1G11870 results in ovule abortion in Arabidopsis [41]. So the different expression of BnaC08g41780D between zws-ms and zws-217 may result some unknown process in fruit development. (3) The knockdown of HAC1 causes disturbances in flower morphology in the Arabidopsis, such as a modified petal shape or even absence of petals [42]. HAC12 has high identity with HAC1, but what morphology changes it can make are still unknown. (4) SEC14-LIKE 3 (ATSFH3) is predominantly expressed in the flowers and involved in the transfer of phosphatidylinositol or phosphatidylcholine phospholipids during flower development in Arabidopsis [43, 44]. In chickpea (Cicer arietinum L.), SEC14-LIKE3 is also considered a potential candidate flowering time-regulating genes by mapping [45]. (5) EC1.3 triggers sperm cell activation during double fertilization and knockdown of it reduces seed number in siliques [46]. (6) The brassinosteroid signal is found involved in ovule initiation and development, by regulating the expression level of BZR1, of which dephosphorylation increases ovule and seed number in Arabidopsis [47]. In this study, BnaA09g44210D has non-synonymous mutation between zws-ms and zws-217, probably causes multi-silique trait by some undiscovered pathways. (7) ATGLK1 (GPRI1) can regulate the silique number in Arabidopsis: the glk1/glk2 double mutant realized only about 1/3 of the mean silique number per plant found in the wild type, due to changes in silique wall and leaf photosynthesis [48]. (8) The Arabidopsis mutant til1–4/til1–4 has altered floral phyllotaxis, reduced ovule number, abnormally developed ovules, and reduced fertility [49]. (9) APA1: In Cynara cardunculus, the aspartic proteinase cardosin B is found playing an important role in ovule function [50]. (10) RLK7: Stamen development depends on a teamwork of receptor-like kinases [51, 52]. (11) CPN10 is involved in response to many physiological process, such as seed abortion in plants [53]. Although no report directly leads to multi-silique trait so far, these orthologs are all related to flower development, flower morphology change, fruit development and even the number of ovules/seeds. It is due to the scarcity of the multi-silique trait to some extent. Moreover, this multi-silique trait is theoretically more complex, since it is controlled by multi-loci. Anyhow, these clues validated the 12 candidate genes.

In summary, the unique multi-silique rapeseed zws-ms, examined in this study represents an enriched rapeseed germplasm resource, which is considered the basis of crop breeding. Thus far, 12 candidate genes in two genomic regions associated with the multi-silique trait have been identified, including seven genes with sequence changes and five genes with expression level differences between the multi-silique and single-silique NILs. This information lays the foundation for future research, such as candidate genes isolation, their functional verification and so on. Several unanswered questions remain. A previous study suggested that the multi-silique is likely controlled by three allelic genes [19], but in the current study, we only identified two associated regions, which harbor no more than two independent genes. This finding suggests that one more gene/factor remains unidentified. We also noted that the multi-silique trait is highly sensitive to the environment. For example, when zws-ms plants were planted in Xindu (altitude of 472 m, 30°47′10″N 104°12′12″E), the multi-silique trait was continuously stable for years; but when they were grown in Ma’erkang (altitude of 2540 m, 31°94′15″N 102°14′29″E), the multi-silique trait disappeared, and all plants had normal siliques (Additional file 12: Fig. S2). Interestingly, once we transferred zws-ms seeds harvested in Ma’erkang back to Xindu, the multi-silique trait was recovered. By contrast, zws-217 maintained the single-silique phenotype at both locations.

Thus, there might be some epigenetic reason for this trait, such as methylation of a base that does not alter the DNA sequence of this underlying locus but is strongly affected by environmental factors. We noticed that there were several genes in the environmental response category in the associated regions with SNP/InDel variations between the two DNA pools: BnaA09g45600D was annotated as “cellular response to water deprivation (GO: 0042631)” and BnaA09g44900D as “response to cold (GO: 0009409)”. Moreover, several DEGs were also annotated to environment response, like BnaA09g45320D was also annotated as “response to cold” and so on. Since these genes are responsive to environmental factors, such as cold, water, and light intensity, their functions vary when their sequences or expression levels change. Thus, they may have different effects under different environmental conditions. Furthermore, Xindu is located in the Sichuan Basin, which has a humid subtropical monsoon climate (with an annual average temperature of 16.2 °C), whereas Ma’erkang is located in a mountainous area in western Sichuan, where the annual average temperature is only 8.6 °C. The environmental factor with the greatest difference between the two locations is temperature. Thus, we propose that temperature is the most important environmental factor affecting the multi-silique trait in zws-ms, although other factors such as light intensity are also likely involved in regulating this trait.

Conclusions

A novel trait in rapeseed, multi-siliques, was investigated at the genome level by association analysis based on WGR, followed by analysis at the transcriptome level. The two regions associated with this trait contain seven genes with non-synonymous mutations or frame-shift mutations annotated to floral organ-related GO terms, as well as five other vital DEGs. These genes are interesting candidate genes for the multi-silique trait.

Methods

Plant materials and population development

The original multi-silique rapeseed line was discovered at the Crop Research Institute, Sichuan Academy of Agricultural Sciences from a population derived from a cross between B. napus and B. rapa. Multi-silique plants were grown and self-pollinated for six successive generations to obtain the homozygous multi-silique zws-ms plants. Single-silique offspring were continuously backcrossed for six generations with zws-ms as the recurrent parent, followed by six continuous generations of self-pollination to obtain the near-isogenic line zws-217 with the single-silique phenotype (Fig. 1). The NILs zws-217 and zws-ms were grown in the field under standard conditions at the Sichuan Academy of Agricultural Sciences in the Xindu District of Sichuan Province, China.

DNA preparation and pool construction

The rate of multi-silique formation per zws-ms plant was determined by phenotyping at the full-bloom stage (BBCH 67). The multi-silique rate per plant was calculated as:

$$ \mathrm{multi}-\mathrm{silique}\ \mathrm{rate}=\frac{\mathrm{number}\ \mathrm{of}\ \mathrm{three}-\mathrm{pistiled}\ \mathrm{flowers}\ \mathrm{per}\ \mathrm{plant}}{\mathrm{total}\ \mathrm{number}\ \mathrm{of}\ \mathrm{flowers}\ \mathrm{in}\ \mathrm{this}\ \mathrm{plant}} $$

Thirty plants with high multi-silique rates were selected for DNA isolation, while 30 zws-217 plants were randomly selected. Fresh leaves were sampled from these plants (zws-ms and zws-217) and subjected to DNA extraction by the cetyltrimethyl ammonium bromide (CTAB) method [54]. RNA contamination was removed from each sample using RNase A. The DNA was quantified using a NanoDrop 2000 Spectrophotometer (Thermo Scientific, USA). Samples with optical density 260/280 ratios ranging from 1.8 to 2.2 were considered adequate. Equal quantities of DNA were pooled by line (multi-silique and single-silique), with 30 samples per line, to a final concentration of 40 ng/μl.

Whole genome re-sequencing

Pooled DNA samples were sheared into ~ 350 bp fragments using a Covaris S2/E210 DNA Shearing kit. The sheared DNA was end-repaired and a single nucleotide (A) overhang was subsequently added to the repaired fragments by adding a Klenow DNA polymerase Fragment (3′ → 5′ exo–) (New England Biolabs, Ipswich, MA, USA) and dATP at 37 °C. Barcodes and Illumina sequencing adapters were ligated to the A-tailed fragments using T4 DNA Ligase (Takara, Dalian, China). PCR was performed on both pooled samples using diluted prepared (sheared and ligated) DNA, deoxyribonucleotide triphosphate, Q5® High-Fidelity DNA Polymerase, and PCR primers. The PCR products were purified using Agencourt AMPure XP Beads (Beckman Coulter, High Wycombe, UK). Fragments 300 to 500 bp in length (including barcodes and adaptors) were excised and purified using a QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany). Gel-purified products were diluted for pair-end sequencing (150 bp on each end) on an Illumina HiSeq system following the standard protocol (Illumina, Inc.; San Diego, CA, USA) at the Biomarker Technologies Corporation (Beijing, China). The sequencing depth was roughly 30 × .

SNP calling

Low-quality reads with quality scores <20e were filtered out and raw reads were sorted according to their barcode sequences. After the barcodes were trimmed, clean HQ reads were mapped to the Brassica_napus.v4.1.fa genome (http://www.genoscope.cns.fr/brassicanapus/data/) using the Burrows-Wheeler Aligner software [55]. SAMtools [56] was used to mark duplicates, and GATK [57] was used for local realignment and base recalibration. An SNP set was formed by combining the results of analysis with GATK and SAMtools via SNP calling using default parameters. SNPs identified between the two lines were regarded as polymorphic for association analysis.

Association analysis

Euclidean distance (ED) association analysis is a method that calculates ED values (quadratic sum root of differences between bulks from the depth of four types of bases) and is satisfied by ED. In theory, the higher the ED value is, the closer the object site [25]. The ED values were calculated as follows [27]:

$$ \mathrm{ED}=\sqrt{{\left({\mathrm{A}}_{\mathrm{ms}}-{\mathrm{A}}_{217}\right)}^2+{\left({\mathrm{C}}_{\mathrm{ms}}-{\mathrm{C}}_{217}\right)}^2+{\left({\mathrm{G}}_{\mathrm{ms}}-{\mathrm{G}}_{217}\right)}^2+{\left({\mathrm{T}}_{\mathrm{ms}}-{\mathrm{T}}_{217}\right)}^2} $$

In this formula, Ams, Cms, Gms, and Tms represent the depth of bases A, C, T, and G on a site in the multi-silique pool, and A217, C217, G217, and T217 represent the depth of bases A, C, T, and G on a site in the single-silique pool, respectively. To remove background noise, the original ED5 value was used, and the adjusted values were fit using the DISTANCE method. Regions over a threshold value were considered to be trait-related candidate regions. The associated threshold value was determined based on ED + 3SD (standard deviation) [27]. InDel-associated regions were obtained via a similar method.

RNA isolation for transcriptome analysis

Three individual plants of the zws-217 (T04, T05, and T06) line were randomly selected for RNA isolation. Eight to ten buds (BBCH 57) per plant were randomly selected and sampled. For the zws-ms, three plants (T01, T02, and T03) were selected, and buds were sampled from each plant. All buds were quick-frozen and stored in liquid nitrogen. Before RNA isolation, the buds from three zws-ms plants were sliced with tweezers and observed. Only buds containing three pistils were considered to be correct multi-silique buds and used for RNA isolation. Buds from each plant were mixed to isolate the RNA, yielding three samples from zws-217 (T04, T05, and T06) and three samples from zws-ms (T01, T02, and T03). Since zws-217 and zws-ms were NILs, and these plants differ from each other only in the multi-silique trait, we reasoned that DEGs between these lines might be related to this trait. Total RNA was isolated using an RNA Isolation Kit (Tiangen, Beijing, China) and the concentration measured using a NanoDrop 2000 (Thermo). RNA integrity was assessed using an RNA Nano 6000 Assay Kit and the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

RNA library preparation and transcriptome sequencing

For each tissue sample, 1 μg RNA was used as input material for RNA sample preparation. Sequencing libraries were generated using a NEBNext Ultra™ RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s instructions with index codes added to associate sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. The remaining overhangs were converted into blunt ends using exonuclease/polymerase. After adenylation of the 3′ ends of DNA fragments, a NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. To select 240 bp cDNA fragments, the library fragments were purified using the AMPure XP system (Beckman Coulter, Danvers, MA USA). The size-selected, adaptor-ligated cDNA was incubated with 3 μl USER Enzyme (NEB, USA) at 37 °C for 15 min, followed by 5 min at 95 °C. PCR was then performed with Phusion High-Fidelity DNA polymerase, universal PCR primers, and an Index (X) Primer. The PCR products were purified (AMPure XP system) and library quality assessed on the Agilent Bioanalyzer 2100 system.

Clustering of index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v4-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced on the Illumina HiSeq X-ten platform and paired-end reads were generated.

Raw reads in fastq format were initially processed using in-house Perl scripts. During this step, clean reads were obtained by removing adapter sequences, reads containing poly-N, and low-quality reads. At the same time, the Q20, Q30, GC-content, and sequence duplication levels of the clean reads were calculated. All downstream analyses were based on clean HQ reads.

Clean reads were mapped to the reference genome sequence, and only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome (http://www.genoscope.cns.fr/brassicanapus/data/). Tophat2 tools (parameters: --read-mismatches 2 --read-edit-dist 2 --max-intron-length 5,000,000 --library-type fr-unstranded --mate-inner-dist 40) were used to map clean reads to the reference genome.

Differential expression analysis

Differential expression analysis of the two lines was performed using the DESeq R package [58]. DESeq provides statistical models based on the negative binomial distribution to identify DEGs using gene expression data. The resulting P values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted fold change of 4 and a p-value ≤0.01 were considered differentially expressed.

Reverse-transcription quantitative PCR (qPCR) validation

To confirm the validity of the DEGs identified by transcriptome sequencing, qPCR was performed. Reverse transcription to cDNA was conducted using a PrimeScript First-strand cDNA Synthesis Kit (Takara, Dalian, China). The relative transcript levels of 10 selected DEGs (at the budding stage BBCH 57) were re-examined by qPCR, with BnACTIN used as the internal control. The primer sequences (Additional file 13: Table S11) for the genes were designed with Premier Primer 3.0 software. All qPCR was performed on a Roche480 instrument (Roche, Basel, Switzerland). Amplification reactions were performed as follows: an initial denaturation step at 95 °C for 3 min, 39 cycles at 95 °C for 10 s, 58 °C for 30 s, and 72 °C for 30 s, a final extension at 72 °C for 10 min and hold at 15 °C.

Abbreviations

BSA:

Bulked-segregant analysis

CTAB:

Cetyltrimethyl ammonium bromide

DEG:

Differentially expressed gene

DH:

Doubled haploid

ED:

Euclidean Distance

GO:

Gene ontology

HQ:

High quality

InDel:

Insertion/deletion

KEGG:

Kyoto Encyclopedia of Genes and Genomes

NGS:

Next-generation sequencing

NIL:

Near-isogenic line

qPCR:

Reverse-transcription quantitative PCR

QTL:

Quantitative trait locus/loci

RAPD:

Random amplified polymorphic DNA

RIL:

Recombinant inbred line

SLAF-seq:

Specific length amplified fragment sequencing

SNP:

Single nucleotide polymorphism

SRAP:

Sequence related amplified polymorphism

SSR:

Simple sequence repeat; WGR: Whole genome re-sequencing

References

  1. 1.

    Zhou HC, Jin L, Li J, Wang XJ. Altered callose deposition during embryo sac formation of multi-pistil mutant (mp1) in Medicago sativa. Genet Mol Res. 2016;15(2). https://doi.org/10.4238/gmr.15027968.

  2. 2.

    Palmer RG, Sandhu D, Curran K, Bhattacharyya MK. Molecular mapping of 36 soybean male-sterile, female-sterile mutants. Theor Appl Genet. 2008;117(5):711–9.

    CAS  Article  Google Scholar 

  3. 3.

    Palmer RG, Horner HT. Genetics and cytology of a genic male-sterile, female-sterile mutant from a transposon-containing soybean population. J Hered. 2000;91(5):378–83.

    CAS  Article  Google Scholar 

  4. 4.

    Fu WQ, Zhao ZG, Ge XH, Ding L, Li ZY. Anatomy and transcript profiling of gynoecium development in female sterile Brassica napus mediated by one alien chromosome from Orychophragmus violaceus. BMC Genomics. 2014;15:61.

    Article  Google Scholar 

  5. 5.

    Li S, Yang L, Deng Q, Wang S, Wu F, Li P. Phenotypic characterization of a female sterile mutant in Rice. J Integr Plant Biol. 2006;48(3):307–14.

    Article  Google Scholar 

  6. 6.

    Arthur L, Ozias-Akins P, Hanna WW. Female sterile mutant in pearl millet: evidence for initiation of Apospory. J Hered. 1993;84(2):112–5.

    Article  Google Scholar 

  7. 7.

    Nair RM, Peck DM, Dundas IS, Samac DA, Moore A, Randles JW. Morphological characterisation and genetic analysis of a bi-pistil mutant (bip) in Medicago truncatula Gaertn. Sex Plant Reprod. 2008;21(2):133–41.

    Article  Google Scholar 

  8. 8.

    Yang Z, Peng Z, Wei S, Liao M, Yu Y, Jang Z. Pistillody mutant reveals key insights into stamen and pistil development in wheat (Triticum aestivum L.). BMC Genomics. 2015;16(1):211.

    Article  Google Scholar 

  9. 9.

    Yang Z, Peng Z, Wei S, Yu Y, Cai P. Identification of differentially expressed genes in three-pistil mutation in wheat using annealing control primer system. Gene. 2011;485(2):81–4.

    CAS  Article  Google Scholar 

  10. 10.

    Beppu K, Kataoka I. Artificial shading reduces the occurrence of double pistils in ‘Satohnishiki’ sweet cherry. Sci Hortic-Amsterdam. 2000;83(3–4):241–7.

    Article  Google Scholar 

  11. 11.

    Malyshev S, Korzun V, Voylokov A, Smirnov V, Rner AB. Linkage mapping of mutant loci in rye (Secale cereale L.). Theor Appl Genet. 2001;103(1):70–4.

    CAS  Article  Google Scholar 

  12. 12.

    Pathak GN, Singh BB. Mutation in Brassica campestris Linn. Var. sarson Prain. Curr Sci. 1949;18:253–4.

    Google Scholar 

  13. 13.

    Xiao L, Zhao H, Zhao Z, Du D, Xu L, Yao Y, Zhao Z, Xing X, Shang G, Zhao H. Genetic and physical fine mapping of a multilocular gene Bjln1 in Brassica juncea to a 208-kb region. Mol Breeding. 2013;32(2):373–83.

    CAS  Article  Google Scholar 

  14. 14.

    Xu P, Lv Z, Zhang X, Wang X, Pu Y, Wang H, Yi B, Wen J, Ma C, Tu J, et al. Identification of molecular markers linked to trilocular gene (mc1) in Brassica juncea L. Mol Breeding. 2014;33(2):425–34.

    CAS  Article  Google Scholar 

  15. 15.

    Katiyar RK, Chamola R, Chopra VL. Tetralocular mustard, Brassica juncea: new promising variability through interspecific hybridization. Plant Breed. 1998;117(4):398–9.

    Article  Google Scholar 

  16. 16.

    Fan C, Wu Y, Yang Q, Yang Y, Meng Q, Zhang K, Li J, Wang J, Zhou Y. A novel single-nucleotide mutation in a CLAVATA3 gene homolog controls a Multilocular silique trait in Brassica rapa L. Mol Plant. 2014;7(12):1788–92.

    CAS  Article  Google Scholar 

  17. 17.

    Guan M, Li X. Effect of 12C heavy ion beams irradiation on rapeseed (Brassica napus). Acta Agron Sin. 2006;32(06):878–84. (in Chinese).

    CAS  Google Scholar 

  18. 18.

    Hu BC, Chen FX, Li QS, Zhang ML. The discovery of multiple pistils CMS in B. campestris L. and its potentiality for breeding. Sci Agric Sin. 1994;27(03):90–1. (in Chinese).

    Google Scholar 

  19. 19.

    Jiang LC, Zhang QX, Pu XB, Zhang ZY. Evaluation of genetic and physiological Characterist ics of aggregate-siliqua rapeseed. Southwest China J Agric Sci. 1998;11(04):62–8 (in Chinese).

    Google Scholar 

  20. 20.

    Michelmore RW, Paran I, Kesseli RV. Identification of markers linked to disease-resistance genes by bulked segregant analysis: a rapid method to detect markers in specific genomic regions by using segregating populations. Proc Natl Acad Sci U S A. 1991;88(21):9828–32.

    CAS  Article  Google Scholar 

  21. 21.

    Shoba D, Manivannan N, Vindhiyavarman P, Nigam SN. SSR markers associated for late leaf spot disease resistance by bulked segregant analysis in groundnut (Arachis hypogaea L.). Euphytica. 2012;188(2):265–72.

    CAS  Article  Google Scholar 

  22. 22.

    Castonguay Y, Cloutier J, Bertrand A, Michaud R, Laberge S. SRAP polymorphisms associated with superior freezing tolerance in alfalfa (Medicago sativa spp. sativa). Theor Appl Genet. 2010;120(8):1611–9.

    CAS  Article  Google Scholar 

  23. 23.

    Laporte V, Merdinoglu D, Saumitou-Laprade P, Butterlin G, Vernet P, Cuguen J. Identification and mapping of RAPD and RFLP markers linked to a fertility restorer gene for a new source of cytoplasmic male sterility in Beta vulgaris ssp. maritima. Theor Appl Genet. 1998;96(6):989–96.

    CAS  Article  Google Scholar 

  24. 24.

    Wang H, Cheng H, Wang W, Liu J, Hao M, Mei D, Zhou R, Fu L, Hu Q. Identification of BnaYUCCA6 as a candidate gene for branch angle in Brassica napus by QTL-seq. Sci Rep. 2016;6:38493.

    CAS  Article  Google Scholar 

  25. 25.

    Geng X, Jiang C, Yang J, Wang L, Wu X, Wei W. Rapid identification of candidate genes for seed weight using the SLAF-Seq method in Brassica napus. PLoS One. 2016;11(1):e147580.

    Article  Google Scholar 

  26. 26.

    Yao Y, Li K, Liu H, Duncan RW, Guo S, Xiao L, Du D. Whole-genome re-sequencing and fine mapping of an orange petal color gene (Bnpc1) in spring Brassica napus L. to a 151-kb region. Euphytica. 2017;213(8):165.

    Article  Google Scholar 

  27. 27.

    Hill JT, Demarest BL, Bisgrove BW, Gorsi B, Su YC, Yost HJ. MMAPPR: mutation mapping analysis pipeline for pooled RNA-seq. Genome Res. 2013;23(4):687–97.

    CAS  Article  Google Scholar 

  28. 28.

    Xia C, Chen L, Rong T, Li R, Xiang Y, Wang P, Liu C, Dong X, Liu B, Zhao D, et al. Identification of a new maize inflorescence meristem mutant and association analysis using SLAF-seq method. Euphytica. 2015;202(1):35–44.

    Article  Google Scholar 

  29. 29.

    Chen W, Yao J, Chu L, Yuan Z, Li Y, Zhang Y. Genetic mapping of the nulliplex-branch gene (gb_nb1) in cotton using next-generation sequencing. Theor Appl Genet. 2015;128(3):539–47.

    CAS  Article  Google Scholar 

  30. 30.

    Qin D, Dong J, Xu F, Guo G, Ge S, Xu Q, Xu Y, Li M. Characterization and fine mapping of a novel barley stage green-revertible albino gene (HvSGRA) by bulked Segregant analysis based on SSR assay and specific length amplified fragment sequencing. BMC Genomics. 2015;16:838.

    Article  Google Scholar 

  31. 31.

    Zheng W, Wang Y, Wang L, Ma Z, Zhao J, Wang P, Zhang L, Liu Z, Lu X. Genetic mapping and molecular marker development for Pi65(t), a novel broad-spectrum resistance gene to rice blast using next-generation sequencing. Theor Appl Genet. 2016;129(5):1035–44.

    CAS  Article  Google Scholar 

  32. 32.

    Xu F, Sun X, Chen Y, Huang Y, Tong C, Bao J. Rapid identification of major QTLs associated with rice grain weight and their utilization. PLoS One. 2015;10(3):e122206.

    Google Scholar 

  33. 33.

    Zhang Z, Li J, Muhammad J, Cai J, Jia F, Shi Y, Gong J, Shang H, Liu A, Chen T, et al. High resolution consensus mapping of quantitative trait loci for fiber strength, length and micronaire on chromosome 25 of the upland cotton (Gossypium hirsutum L.). PLoS One. 2015;10(8):e135430.

    Google Scholar 

  34. 34.

    Varshney RK, Nayak SN, May GD, Jackson SA. Next-generation sequencing technologies and their implications for crop genetics and breeding. Trends Biotechnol. 2009;27(9):522–30.

    CAS  Article  Google Scholar 

  35. 35.

    Hu MJ, Zhang HP, Liu K, Cao JJ, Wang SX, Jiang H, Wu ZY, Lu J, Zhu XF, Xia XC, et al. Cloning and characterization of TaTGW-7A gene associated with grain weight in wheat via SLAF-seq-BSA. Front Plant Sci. 2016;7:1902.

    PubMed  PubMed Central  Google Scholar 

  36. 36.

    Wei QZ, Fu WY, Wang YZ, Qin XD, Wang J, Li J, Lou QF, Chen JF. Rapid identification of fruit length loci in cucumber (Cucumis sativus L.) using next-generation sequencing (NGS)-based QTL analysis. Sci Rep. 2016;6:27496.

    CAS  Article  Google Scholar 

  37. 37.

    Fallath T, Kidd BN, Stiller J, Davoine C, Bjorklund S, Manners JM, Kazan K, Schenk PM. MEDIATOR18 and MEDIATOR20 confer susceptibility to Fusarium oxysporum in Arabidopsis thaliana. PLoS One. 2017;12(4):e176022.

    Article  Google Scholar 

  38. 38.

    Komeda Y. Genetic regulation of time to flower in Arabidopsis Thaliana. Annu Rev Plant Biol. 2004;55(1):521–35.

    CAS  Article  Google Scholar 

  39. 39.

    Mouradov A, Cremer F, Coupland G. Control of flowering time: interacting pathways as a basis for diversity. Plant Cell. 2002;14 Suppl:S111–30.

    CAS  Article  Google Scholar 

  40. 40.

    Chen L, Cheng JC, Castle L, Sung ZR. EMF genes regulate Arabidopsis inflorescence development. Plant Cell. 1997;9(11):2011–24.

    CAS  Article  Google Scholar 

  41. 41.

    Berg M, Rogers R, Muralla R, Meinke D. Requirement of aminoacyl-tRNA synthetases for gametogenesis and embryo development in Arabidopsis. Plant J. 2005;44(5):866–78.

    CAS  Article  Google Scholar 

  42. 42.

    Boycheva I, Vassileva V, Revalska M, Zehirov G, Iantcheva A. Different functions of the histone acetyltransferase HAC1 gene traced in the model species Medicago truncatula, Lotus japonicus and Arabidopsis thaliana. Protoplasma. 2017;254(2):697–711.

    CAS  Article  Google Scholar 

  43. 43.

    Mo P, Zhu Y, Liu X, Zhang A, Yan C, Wang D. Identification of two phosphatidylinositol/phosphatidylcholine transfer protein genes that are predominately transcribed in the flowers of Arabidopsis thaliana. J Plant Physiol. 2007;164(4):478–86.

    CAS  Article  Google Scholar 

  44. 44.

    Wang H, You C, Chang F, Wang Y, Wang L, Qi J, Ma H. Alternative splicing during Arabidopsis flower development results in constitutive and stage-regulated isoforms. Front Genet. 2014;5:25.

    PubMed  PubMed Central  Google Scholar 

  45. 45.

    Upadhyaya HD, Bajaj D, Das S, Saxena MS, Badoni S, Kumar V, Tripathi S, Gowda CLL, Sharma S, Tyagi AK, et al. A genome-scale integrated approach aids in genetic dissection of complex flowering time trait in chickpea. Plant Mol Biol. 2015;89(4):403–20.

    CAS  Article  Google Scholar 

  46. 46.

    Sprunck S, Rademacher S, Vogler F, Gheyselinck J, Grossniklaus U, Dresselhaus T. Egg cell-secreted EC1 triggers sperm cell activation during double fertilization. Science. 2012;338(6110):1093–7.

    CAS  Article  Google Scholar 

  47. 47.

    Huang HY, Jiang WB, Hu YW, Wu P, Zhu JY, Liang WQ, Wang ZY, Lin WH. BR signal influences Arabidopsis ovule and seed number through regulating related genes expression by BZR1. Mol Plant. 2013;6(2):456–69.

    CAS  Article  Google Scholar 

  48. 48.

    Zhu X, Zhang L, Kuang C, Guo Y, Huang C, Deng L, Sun X, Zhan G, Hu Z, Wang H, et al. Important photosynthetic contribution of silique wall to seed yield-related traits in Arabidopsis thaliana. Photosynth Res. 2018;137(3):493–501.

    CAS  Article  Google Scholar 

  49. 49.

    Jenik PD, Jurkuta REJ, Barton MK. Interactions between the cell cycle and embryonic patterning in Arabidopsis uncovered by a mutation in DNA polymerase. Plant Cell Online. 2005;17(12):3362–77.

    CAS  Article  Google Scholar 

  50. 50.

    Figueiredo R, Duarte P, Pereira S, Pissarra J. The embryo sac of Cynara cardunculus: ultrastructure of the development and localisation of the aspartic proteinase cardosin B. Sex Plant Reprod. 2006;19(2):93–101.

    CAS  Article  Google Scholar 

  51. 51.

    Cai W, Zhang D. The role of receptor-like kinases in regulating plant male reproduction. Plant Reprod. 2018;31(1):77–87.

    CAS  Article  Google Scholar 

  52. 52.

    Zhao D. Control of anther cell differentiation: a teamwork of receptor-like kinases. Sex Plant Reprod. 2009;22(4):221–8.

    CAS  Article  Google Scholar 

  53. 53.

    Hanania U, Velcheva M, Or E, Flaishman M, Sahar N, Perl A. Silencing of chaperonin 21, that was differentially expressed in inflorescence of seedless and seeded grapes, promoted seed abortion in tobacco and tomato fruits. Transgenic Res. 2007;16(4):515–25.

    CAS  Article  Google Scholar 

  54. 54.

    Doyle J, Doyle J. Isolation of plant DNA from fresh tissues. Focus. 1990;12:13–5.

    Google Scholar 

  55. 55.

    Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    CAS  Article  Google Scholar 

  56. 56.

    Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  Google Scholar 

  57. 57.

    McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303.

    CAS  Article  Google Scholar 

  58. 58.

    Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–8.

    Article  Google Scholar 

Download references

Acknowledgements

Not Applicable

Funding

This research was supported by: The National Key Research and Development Plan (2016YFD0101305); Modern Agro-industry Technology Research System of China (CARS-12);

Major Science and Technology Special Subject of Sichuan Province (2018NZDZX0003); Scientific Observing and Experimental Station of Oil Crops in the Upper Yangtze River, Ministry of Agriculture, P. R. China (09203020); Financial Innovation Ability Promotion Project of Sichuan Province (2016ZYPZ-013); The National Key Research and Development Plan (2018YFD0100500); National Natural Science Foundation of China (31560402). The funders were not involved in the study design or collection, analysis, or interpretation of the data.

Availability of data and materials

The reference genome sequence was based on the Brassica napus L. reference genome V4.1 (http://www.genoscope.cns.fr/brassicanapus/data/); the additional (supplementary) files are attached; Raw data obtained by WGR and transcriptome sequencing were deposited in the Sequence Read Archive (SRA) database in NCBI under accession number PRJNA491622.

Author information

Affiliations

Authors

Contributions

LC and JZ carried out the experiments, performed data analysis, and drafted the manuscript; KL, LW, and HW revised the manuscript and optimized the figures; HL, JJ, CC, and BZ assisted with DNA/RNA extraction and data analysis; LJ designed and supervised the experiments; All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Liangcai Jiang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not Applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S1. Whole genome with associated regions. (DOCX 209 kb)

Additional file 2:

Table S1. Summary of the transcriptome sequencing data. (DOCX 15 kb)

Additional file 3:

Table S2. Information about the mapped reads based on the transcriptome sequencing data. (DOCX 15 kb)

Additional file 4:

Table S3. Annotations of the 2041 genes in the two intersection associated regions. (XLSX 340 kb)

Additional file 5:

Table S4. GO enrichment in biological processes (BP) category of the genes. (XLSX 123 kb)

Additional file 6:

Table S5. GO enrichment in cellular components (CC) category of the genes. (XLSX 26 kb)

Additional file 7:

Table S6. GO enrichment in molecular function (MF) category of the genes. (XLSX 50 kb)

Additional file 8:

Table S7. GO annotations of the 104 genes with non-synonymous mutation SNPs in the associated regions. (DOCX 34 kb)

Additional file 9:

Table S8. GO annotations of the 19 genes with frame-shift mutation InDels in the associated regions. (DOCX 17 kb)

Additional file 10:

Table S9. KEGG enrichment of the 2041 genes in the intersection associated regions. (XLSX 12 kb)

Additional file 11:

Table S10. Information about the DEGs between zws-ms and zws-217. (DOCX 40 kb)

Additional file 12:

Figure S2. Stability of the multi-silique trait under one environment and variation of this trait across different environments. (DOCX 1840 kb)

Additional file 13:

Table S11. Primers used for qPCR. (DOCX 15 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Chai, L., Zhang, J., Lu, K. et al. Identification of genomic regions associated with multi-silique trait in Brassica napus. BMC Genomics 20, 304 (2019). https://doi.org/10.1186/s12864-019-5675-4

Download citation

Keywords

  • Association analysis
  • Brassica napus L.
  • Multi-silique
  • Near-isogenic line
  • Transcriptome sequencing
  • Whole genome re-sequencing