- Research Article
- Open Access
An integrated analysis of QTL mapping and RNA sequencing provides further insights and promising candidates for pod number variation in rapeseed (Brassica napus L.)
BMC Genomics volume 18, Article number: 71 (2017)
As the most important yield component in rapeseed (Brassica napus L.), pod number is determined by a series of successive growth and development processes. Pod number shows extensive variation in rapeseed natural germplasm, which is valuable for genetic improvement. However, the genetic and especially the molecular mechanism for this kind of variation are poorly understood. In this study, we conducted QTL mapping and RNA sequencing, respectively, using the BnaZNRIL population and its two parental cultivars Zhongshuang11 and No.73290 which showed significant difference in pod number, primarily due to the difference in floral organ number.
A total of eight QTLs for pod number were identified using BnaZNRIL population with a high-density SNP linkage map, each was distributed on seven linkage groups and explained 5.8–11.9% of phenotypic variance. Then, they were integrated with those previously detected in BnaZNF2 population (deriving from same parents) and resulted in 15 consensus-QTLs. Of which, seven QTLs were identical to other studies, whereas the other eight should be novel.
RNA sequencing of the shoot apical meristem (SAM) at the formation stage of floral bud primordia identified 9135 genes that were differentially expressed between the two parents. Gene ontology (GO) analysis showed that the top two enriched groups were S-assimilation, providing an essential nutrient for the synthesis of diverse metabolites, and polyamine metabolism, serving as second messengers that play an essential role in flowering genes initiation. KEGG analysis showed that the top three overrepresented pathways were carbohydrate (707 genes), amino acid (390 genes) and lipid metabolisms (322 genes).
In silico mapping showed that 647 DEGs were located within the confidence intervals of 15 consensus QTLs. Based on annotations of Arabidopsis homologs corresponding to DEGs, nine genes related to meristem growth and development were considered as promising candidates for six QTLs.
In this study, we discovered the first repeatable major QTL for pod number in rapeseed. In addition, RNA sequencing was performed for SAM in rapeseed, which provides new insights into the determination of floral organ number. Furthermore, the integration of DEGs and QTLs identified promising candidates for further gene cloning and mechanism study.
Pod number is one of the three yield components (pod number per plant, seed number per pod and seed weight) and breeding targets in rapeseed (Brassica napus L.). Among the three components of yield in rapeseed, pod number shows the highest correlation with yield, which suggests it to be the major contributor to yield . Pod number, especially from branch inflorescence and whole plant, is also the most variable among the three yield components, in accordance with its modest heritability observed in most studies [2, 3]. In rapeseed germplasm resources, pod number shows extensive natural variation, which is invaluable for the genetic improvement . However, the genetic and especially the molecular mechanism for this kind of variation are poorly understood.
Pod number is also a very complex trait that is multiplicatively determined by its three components: the number of flowers differentiated, the proportion of ovaries successfully fertilized, and the rate of fertilized ovaries developed into pods, which are determined by the flower bud differentiation, fertilization and pod development, respectively. Of these, flower bud differentiation is the first critical developmental stage and morphogenesis process which determines the number of floral organs . This process is undertaken with the interaction of internal (such as carbohydrates, phytohormones, polyamine etc.) and external (such as light, temperature, water, and fertilizer etc.) factors . Whereas, genetic factors are the most important control combined with internal signals (also genetically controlled) and influenced by environmental interactions. During the reproductive phase, the shoot apical meristem (SAM) produces inflorescence meristem (IM) that quickly develops into floral meristems (FMs) that, in turn, produce floral primordia . In Arabidopsis, more than 100 regulators had been characterized to involve in SAM growth and development [6–10]. In addition, a complex regulatory network of hormones, transcription factors, enzymes, microRNA and other cellular components had been developed . In rice, several genes involving in the regulation of floral organ number had been isolated, such as FON1-FON4 [12–15]. However, molecular mechanism of floral bud differentiation remained relatively unclear in Brassica napus.
Pod number is a typical quantitative trait, which shows continuous variation and is very sensitive to the environmental conditions . To the present, more than 80 pod number quantitative trait loci (QTLs) have been identified from nearly ten linkage mapping populations . However, only a few of these QTLs were repeatedly detected, in accordance with the modest heritability of pod number. In addition, almost all of these pod number QTLs showed a moderate effect . Therefore, it is difficult to narrow down these pod number QTLs and identify the underlying candidate genes.
The previous expression profiling studies largely relied on hybridization-based technologies, which was unable to fully catalogue and quantify the diverse RNA molecules that are expressed from genomes over a wide range of levels . With the rapid advancement and decreased price of high-throughput sequencing technology in recent years, RNA sequencing has been widely applied in various species to conduct annotation and quantification of all genes and their isoforms across samples . However, a large amount of differentially expressed genes (DEGs) were usually identified between contrast samples through RNA sequencing. Therefore, integration of DEGs and QTLs is considered to be a promising method to identify potential candidate genes .
In our previous study, we screened two rapeseed cultivars Zhongshuang11 and No.73290 that showed significant difference in pod number  mainly due to the difference in floral organ number. The main objectives of our study are to: (1) dissect the genetic mechanism of pod number variation by QTL mapping using BnaZNRIL and BnaZNF2 population derived from Zhongshuang11 and No.73290; (2) dissect the molecular mechanism of pod number variation by characterizing the expression profile and DEGs of SAM for Zhongshuang11 and No.73290 using RNA sequencing technology; (3) integrate the QTLs and DEGs to identify potential candidate genes for further functional and mechanism studies.
Mapping of QTLs for pod number using BnaZNRIL population
The pod number of No.73290 was much larger than that of Zhongshuang11 in all of the three investigated environments (Table 1). The pod number of the BnaZNRIL population showed normal or near-normal distribution in the four environments (Fig. 1), indicating a quantitative inheritance suitable for QTL mapping . The broad-sense heritability was 0.67, 0.62 and 0.62 for pod number of main inflorescence, branch inflorescence and whole plant, respectively (Additional file 1: Table S1).
A high-density linkage map comprising 2264 SNP markers was constructed using the BnaZNRIL population. After deleting two non-reproducible suggestive QTL, a total of eight QTLs were identified for pod number in four experiments, six and two for main inflorescence (PNm) and branch inflorescence (PNb), respectively (Fig. 2). The identified QTLs were distributed on A01 (qPN.A01-2, −3,), A02 (qPN.A02-1), A03 (qPN.A03-3), A04 (qPN.A04-1), A06 (cqPN.A06-1), C04 (qPN.C04-1) and C06 (qPN.C06-2) linkage groups. These QTLs for PNb and PNm explained 9.3–9.9% and 5.8–11.9% of phenotypic variance, respectively. The additive effect of these QTLs for PNb and PNm ranged from −16.3 to −13.8 and from −4.1 to 2.5, respectively. Furthermore, the corresponding genomic intervals of the eight QTLs were determined through mapping the probe sequences of SNP markers within their confidence intervals to the Brassica napus genome (Additional file 2: Table S2).
Because QTL mapping of pod number had been previously performed using the BnaZNF2 population derived from the same parents , meta-analysis was carried out to integrate pod number QTLs detected in the two studies. After meta-analysis, a total of 15 consensus-QTLs were obtained (Table 2), of which seven were identical to the previously identified ones, whereas the other eight should be novel (Additional file 3: Table S3).
Because the flowering time of the two parents (Zhongshuang11 and No.73290) for both BnaZNF2 and BnaZNRIL populations differed for about one week, therefore the 15 pod number consensus-QTLs were also compared with the flowering time QTLs detected in both populations. The results showed that most of the QTLs for both traits were not overlapped (Shi et al. unpublished data).
Based on rapeseed genome annotation , a total of 6164 genes resided in regions of 15 QTLs. Number of rapeseed genes in QTL region ranged from 63 (qPN.A02-1) to 1103 (qPN.C04-1), with the mean of 411 (Table 3). The majority of these genes have never been functionally annotated in Brassica napus.
Transcriptome sequencing and mapping
To dissect the possible molecular mechanism of pod number difference between Zhongshuang11 and No.73290, transcriptome sequencing was performed using the shoot apical meristem (SAM) of Zhongshuang11 and No.73290 because their pod number difference was mainly attributable to the floral organ number. Total RNA, isolated from SAM of ten individuals of Zhongshuang11 and No.73290 at the formation stage of floral bud primordia, was pooled and transcribed into cDNA. Then, the two cDNA libraries were constructed and sequenced via Illumina Hiseq 2000 platform. Consequently, 172,006,032 and 170,055,032 raw reads were generated for Zhongshuang11 and No.73290, respectively. After removal of sequences with adaptor, duplication and low quality, we obtained 157,471,012 and 155,457,898 high-quality clean reads for Zhongshuang11 and No.73290, respectively. A summary of the transcriptome sequencing is exhibited in Table 4.
A total of 127,275,824 (80.82%), and 125,875,780 (80.97%) clean reads for Zhongshuang11 and No.73290, respectively, were successfully mapped to the released reference genome  of Brassica napus using SOAPaligner/soap2  (Table 5). Of which, 106,578,237 (67.68%) and 105,439,605 (67.83%) clean reads for Zhongshuang11 and No.73290, respectively, were uniquely mapped to the reference genome.
Identification and validation of DEGs between Zhongshuang11 and No.73290
All the uniquely mapped reads were selected for further analysis. To remove technical biases inherent in the sequencing process, most notably the length of the RNA species and the sequencing depth of samples, RPKM (reads per kilobase per million reads) was used to normalize these measures . A total of 68030 (67.33%) and 69733 (69.02%) genes expression were detected for Zhongshuang11 and No.73290 respectively. Notably, 64440 (87.89%) genes expression was commonly detected in both Zhongshuang11 and No.73290 (Fig. 3a). Genes with average RPKMs beyond 60 were selected to conduct gene ontology (GO) analysis (Fig. 3b). These genes involved in biological processes including electron transport or energy pathway, protein metabolism and other biological process. The most enriched category in molecular function was structural molecule activity. The most overrepresented category in cellular component were ribosome, followed by cytosol, other cellular components and cell wall.
Using false discovery rate (FDR) < 0.001 and fold change > 2 as the significant level, a total of 9135 differentially expressed genes (DEGs) were identified. Of which, 5507 (60.28%) and 3628 (39.72%) DEGs were respectively up and down regulated, in No.73290 compared to Zhongshuang11 (Additional file 4: Table S4).
To validate the results of RNA sequencing, qRT-PCR was conducted on 12 randomly selected DEGs with 7 down-regulated and 5 up-regulated genes. As shown in Fig. 4, the expression patterns for 12 DEGs were generally consistent between RNA sequencing and qRT-PCR data. Whereas, differences in fold changes measured by RNA sequencing and qRT-PCR were observed, which is consistent with other studies .
Functional annotation of the DEGs
All the DEGs were annotated based on their similarity (E value ≤1e-5) in six public databases, including GO, KEGG, Pfam, InterPro, SwissProt, and TrEMBLE (Table 6). In total, 8479 DEGs (92.82%) have significant similarity in at least one database.
Because the majority of the rapeseed reference genes have no functional annotation. These DEGs were then annotated by sequence alignment with those in TAIR 10 (E-value cutoff of 1e-5). Finally, 8305 (90.91%) of 9135 DEGs matched at least one gene in Arabidopsis (Additional file 5: Table S5). To overview functions of DEGs, the 8305 annotated DEGs were assigned to at least one Gene Ontology (GO) category that belonged to three major terms: cellular component, molecular function and biological process. After the absolute gene numbers in each group were normalized to the frequency of group over all Arabidopsis genes (Fig. 5), S-assimilation was the most over-represented group, indicating that the role of S-assimilation in SAM growth and development. In plant, the uptake of sulfate and its assimilation provides an essential nutrient for the synthesis of diverse metabolites, including cystesine, methionine, glutathione, vitamin cofactors and so on . These metabolites affect carbon nitrogen ratio, which likely contributes to floral bud development . The following over-represented group was polyamine metabolism. Previous studies indicated that polyamine served as second messengers playing an essential role in flowering genes initiation . Other significantly enriched groups, such as N-metabolism, amino acid metabolism, lipid metabolism or hormone metabolism, displayed the active metabolic status of SAM.
To further study the biological pathways that related to SAM growth and development, the pathway enrichment analysis was conducted in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. These DEGs were mapped to 190 KEGG pathways (Additional file 6: Table S6), which were grouped into five categories: cellular processes, genetic information processing, environmental information processing, metabolism, and organismal systems (Fig. 6). Of which, metabolism (2835 genes, 72.01%) was most enriched, which indicated that developmental differences of SAM between Zhongshuang11 and No.73290 were largely related to metabolism. Moreover, the top three represented pathways were carbohydrate metabolism (707 genes), amino acid metabolism (390 genes) and lipid metabolism (322 genes), which all belonged to the metabolism group. This is understandable, because the flower bud differentiation is dependent first on the nutrient level in the body of plant, which is reflected by the cytosol concentration (in the shoot apex growing point) that is determined by the metabolic process. Carbohydrate (as the structure and energy matter) accumulation is closely related to flower bud differentiation . Increasing amino acid content promotes flower bud formation . Lipid metabolism contributes to cell membrane and other parts. Because the SAM is a reservoir of undifferentiated stem cells that function as a continuous source of new cells . These results provide further insight into the molecular mechanism responsible for floral organ number variation in rapeseed.
Integration of DEGs with QTLs
To further understand the roles of these DEGs played in regulating SAM development and the final flower/pod number, they were integrated with the total of 15 pod number consensus-QTLs identified in both BnaZNRIL and BnaZNF2 population by in silico mapping. A total of 647 DEGs were located in the regions of the 15 QTLs, and the DEGs number within each QTL ranged from 11 (qPN.A02-1) to 74 (qPN.A01-2/ qPN.A01-3), with a mean of 43 (Table 3). Of which, 604 DEGs (93.35%) showed significant homology with Arabidopsis genes and the most have functional annotation.
Interestingly, nine DEGs underlying six QTLs were known to play an important role in regulating meristem growth and development. The nine DEGs resided in the regions of qPN.A01-3(BnaA01g22100D), qPN.A03-1(BnaA03g25890D), qPN.A03-2(BnaA03g29180D, BnaA03g29810D), qPN.A05-1(BnaA05g12220D), qPN.C02-1(BnaC02g02900D, BnaC02g03640D), and qPN.C06-2(BnaC06g28880D, BnaC06g29980D) (Table 7). Their functions included transcription factor, auxin receptor, enzyme, and other cellular elements. For example, BnaA03g25890D (underlying qPN.A03-1) is homologous to Arabidopsis ABP1 (AT4G02980), which is an extra-cellar auxin receptor that involves in the maintenance of anisotropic growth to initiate new floral primordial . BnaC06g29980D (underlying qPN.C06-2) is homologous to Arabidopsis AP1 (AT1G69120), a floral homeotic gene encoding MADS-domain transcription factor, which plays an essential role in floral meristem determinacy, maintenance of floral meristem identity, flower development, and the transcriptional activation of several flowering time genes including AG, SVP, SOC1 and AGL24 [29, 30]. BnaC02g02900D (underlying qPN.C02-1) is homologous to Arabidopsis TFL1 (AT5G03840), which is involved in floral initiation process and control floral meristem identity . The nine genes should be considered as promising candidates that could be used for further study. Our study also indicated that integration of transcriptome sequencing with QTL mapping could be an efficient method to narrow down the number of candidate genes within the QTL region.
Novel QTLs identified for pod number
In the present study, a total of eight QTLs for pod number were identified in the BnaZNRIL population, and were integrated with those previously detected in BnaZNF2 population , which resulted in a total of 15 consensus-QTLs. These QTLs were distributed on ten linkage groups (A01, A02, A03, A04, A05, A06, A09, C02, C04 and C06). Most of these QTLs showed a moderate effect (R 2 < 10%) and only one (cqPN.A06-1) can be considered as a “major” QTL. Of these, only four QTLs were repeatedly detected and the other 11 were specific (Table 2). This strongly indicated that pod number is very sensitive to environmental condition, and accorded well with the moderate heritability of pod number found in the current and previous studies [2, 3].
More than 70 QTLs for pod number had been detected in the other studies, which were distributed on 17 (excluding C04 and C07) of the 19 linkage groups . To accurately identify the positional relationship of pod number QTLs detected in our and other studies, comparative QTL analysis was performed based on the physical map of Brassica napus . Of the total of 89 pod number QTLs identified in our and other studies, 83 ones were successfully mapped on the physical map. This represents the first relatively comprehensive genetic architecture of pod number in rapeseed. Of the 15 consensus-QTLs detected in our studies, the genomic regions of seven QTLs (qPN.A01-2, qPN.A02-1, qPN.A03-1, qPN.A03-2, qPN.A03-3, qPN.C06-1 and qPN.C06-2) were overlapped with those reported in other studies, whereas the other eight (qPN.A01-3, cqPN.A01-1, qPN.A04-1, qPN.A05-1, cqPN.A06-1, qPN.A09-1, qPN.C02-1 and qPN.C04-1) should be novel. More importantly, we detected a pod number QTL (qPN.C04-1) on C04 linkage group for the first time.
Characteristics of SAM DEGs
Transcriptome sequencing of SAM had been documented in several crops such as maize  and soybean , but had rarely been reported in Brassica napus. To our knowledge, this is the second report on SAM transcriptome sequencing in rapeseed. In this study, KEGG analysis showed that the most enriched group of these DEGs belong to metabolism, especially carbohydrate metabolism, which provides the energy and material basis to produce the difference of the number of differentiated flower buds and opening flowers. We proposed that these carbohydrates should be mostly transferred from the leaf, because it is the major source of photosynthesis during the vegetative and early reproductive stage.
To date, a total of 2296 transcription factors have been identified and classified into 58 families in Arabidopsis (http://planttfdb.cbi.pku.edu.cn/index.php?sp=Ath). Based on the functional annotation of corresponding Arabidopsis homologous of these DEGs, nearly all transcription factors families (52 of the total 58 families) were found. This indicated that transcription factors may play an important role in SAM development. In addition, phytohormones metabolism was overrepresented. The roles of several types of hormones in SAM development have been extensively studied. Of which, cytokinins and auxins act as two major hormones to involve in meristem function and maintenance . Auxins are a positive regulator to induce lateral organ initiation, such as flower bud, leaf. Cytokinins play a role in meristem maintenance and in controlling meristematic properties, such as cell proliferation. Furthermore, auxins interact with cytokinins to regulate SAM development [34, 35]. In the SAM, gibberellin activity is down-regulated and its activity is antagonistic with cytokinins activity . In the SAM, brassinosteroids play a role in the spatiotemporal control of organ boundary formation and morphogenesis . What’s more, phytohormones and transcription factors can cooperate to balance meristem maintenance and organ production . These results showed a complex regulatory network involved in SAM development.
Candidates for pod number
Previous studies into pod number mainly focused on QTL mapping (in genomics level) and no study has been conducted in transcriptomics. RNA sequencing is a new high-throughput, high-sensitivity and high-speed technology, which is widely used for transcriptome profiling . In comparison with the previous methods, RNA sequencing possesses three key advantages: first, RNA sequencing is not limited to detect transcripts that correspond to existing genomic sequences; Second, RNA sequencing has very low background signal; Third, transcripts detected by RNA sequencing have a large dynamic range of expression levels . The integration of the transcriptome profiling and QTL mapping is a very useful and popular strategy toward discovery of candidate genes .
Our previous study showed that pod number of No.73290 was about half greater than that of Zhongshuang11 , which was mainly due to their difference in floral organ number according to results of successive observation in several years. According to the relevant research, floral organ number was mainly determined by the SAM development [12, 39]. Thus, SAM of both Zhongshuang11 and No.73290 was characterized using the RNA sequencing technology and the identified DEGs was combined with QTLs to narrow the candidates.
According to the functional annotation of the corresponding Arabidopsis homologues of these DEGs, nine involved in meristem growth and development were considered to be the candidates, which should be important targets for further functional validation.
Transcription factors (such as WUS, STM) play a central role in SAM growth and development [11, 40]. Of the above-mentioned nine candidates, the Arabidopsis homologues (AP1, BLH9, and AT1G60340) of BnaC06g29980D, BnaC02g03640D and BnaA01g22100D were transcription factor-encoding genes. AP1, a MAD box transcription factor, activates floral organ identity genes to promote FMs formation by interplaying with LFY [29, 41, 42]. Loss function of AP1 prevents the formation of flowers in the axils of sepal to form FMs . In our study, the expression level of BnaC06g29980D was higher in No.73290 than in Zhongshuang11, which might promote flower formation. BLH9, a homeodomain transcription factor, contributes to spatial expression patterns of boundary genes and floral induction by FT . Loss function of this gene impairs stem cell maintenance and blocks internode elongation and flowering . AT1G60340 belongs to the NAC family, whose proteins have a consensus sequence known as the NAC domain (petunia NAM and Arabidopsis ATAF1, ATAF2, and CUC2) . Previous studies showed that mutated NAM (NO APICAL MERISTEM) genes failed to form SAM in petunia and that mutated CUC2 causes defect in the formation of the SAM in Arabidopsis .
Phytohormones, such as auxin, cytokinins and gibberellin etc., play an important role in the shoot apical and floral meristem functions . Of the above-mentioned nine candidates, BnaA03g25890D is homologous to Arabidopsis ABP1, which functions as an auxin receptor and its knockdown leads to an enhanced degradation of AUX/IAA repressors . In current study, the expression level of BnaA03g25890D was lower in No.73290 than in Zhongshuang11, which might promote flower primordium initiation.
BnaA03g29180D is homologous to ATSK12, which encodes SHAGGY-like protein kinase involved in meristem organization . BnaA05g12220D is homologous to ASL5, which encodes the protein containing a conserved LOB domain. The T-DNA insertion mutant of ASL5 displayed lost apical dominance, abnormal inflorescence compare to wide type [49, 50]. In our study, the expression level of BnaA03g29180D and BnaA05g12220D in No.73290 was higher more than in Zhonghsuang11, which suggested that BnaA03g29180D and BnaA05g12220D may positively regulate the SAM development.
BnaC02g02900D is homologous to TFL1 (encoding a small transcription cofactors), which controls inflorescence meristem identity and is involved in the floral initiation process [31, 51]. To date, the function of TFL1 has been extensively studied in soybean , rose  and black cherry  etc., which all indicated its role in the transcription repression in floral development. In our study, the expression level of BnaC02g02900D was lower in No.73290 than in Zhongshuang11, which might promote floral development. BnaA03g29810D is homologous to NSN1, which encodes a nucleolar GTP-binding protein and is required for maintenance of inflorescence meristem identity and floral organ development.. Loss function mutant of NSN1 showed a smaller inflorescence meristem dome in comparison to the wide type plants at young stages . In current study, BnaA03g29810D showed higher expression level in No.73290 than in Zhongshuang11, which might promote floral organ development. BnaC06g28880D is homologous to TEL2, which is similar to the terminal ear1 (TE1) gene in maize. The te1 mutant displayed earlike appearance of the terminal inflorescence, which indicated that TEL2 may contribute to normal inflorescence and lead to floral development . The expression level of BnaC06g28880D was higher in No.73290 than in Zhongshuang11, which might promote inflorescence and flower development.
To further validate the nine candidate genes, their over-expression and knockout vector have been constructed and transformed in our laboratory.
To investigate the genetic and molecular mechanism of pod number variation in rapeseed, we conducted QTL mapping and RNA sequencing, respectively, using the BnaZNRIL/BnaZNF2 populations and its two parents Zhongshuang11 and No.73290. As a result, 15 consensus-QTLs were identified through meta-analysis. Go and KEGG analysis of 9135 DEGs between the SAM of Zhongshuang11 and No.73290 provided new and further insights into pod number variation. An integrated analysis of QTLs and DEGs identified several promising candidates for these QTLs, which laid a solid foundation for further functional and mechanism research.
Population construction, field experiments and trait investigation
Two linkage populations, BnaZNF2 and BnaZNRIL, were used for mapping and integration of QTLs for pod number in our study. Both populations were derived from the two sequenced rapeseed cultivars, Zhongshuang11 (de novo sequencing) and No.73290 (re-sequencing) that showed significant difference on pod number. The BnaZNF2 population included 184 F2 individuals, which had been reported previously . The BnaZNRIL population included 184 F7 generation RILs (recombinant inbred lines) derived from the above-mentioned 184 F2 individuals using single-seed descent.
The BnaZNRIL population was planted in Wuhan (Hu Bei province, China) and Zhengzhou (He Nan province, China) from Oct. 2012 to May 2013 and from Oct. 2013 to May 2014 (code W13RIL, Z13RIL, W14RIL, and Z13RIL, respectively) followed a randomized complete block design with two replications. Each block contained three rows, with 33 cm distance between rows and 16.7 cm distance between individuals. The field management followed standard agriculture practice. At maturity, 10 representative plants from the middle of the second row of each block were harvested.
Effective pod number was investigated according to the previous study . Pod number included three categories: pods from the main inflorescence (PNm), branch inflorescence (PNb) and whole plant (PNw), respectively.
SNP genotyping and linkage map construction
To genotype the BnaZNRIL population, The Brassica 60 K Illumina® Infinium SNP array was recently developed by the international Brassica Illumina SNP consortium. The array hybridization and data processing were carried out in Emei Tongde Co. (Beijing) according to the manufacturer’s recommendations. Leaf tissue from seedlings of BnaZNRIL population was used to extract genomic DNA according to the CTAB method .
The genetic linkage map was constructed using the software JoinMap 4.1 (https://www.kyazma.nl/index.php/JoinMap/) with a threshold for goodness-of-fit ≤5, recombination frequency of < 0.4 and minimum logarithm of odds (LOD) score of 2.0. The genetic distances were measured based on the Kosambi function. To avoid the potential errors, double-crossover events were checked.
QTL mapping and meta-analysis
QTL mapping was performed using the composite interval mapping (CIM) method incorporated into WinQTLCart v2.5 software (http://statgen.ncsu.edu/qtlcart/WQTLCart.htm). The parameters including walk speed, number of control markers, window size and regression method were set to 1 cM, 5, 10 cM, and forward regression, respectively. Permutation analysis  with 1000 repetitions was used to determine the LOD threshold. The LOD value corresponding to P = 0.05 (3.0–4.6) was used to detected significant QTLs. Moreover, a lower LOD value corresponding to P = 0.10 (2.6–4.1) was employed to identify suggestive QTLs with small effects.
Meta-analysis was used to estimate the number and position of the true QTLs underlying the QTLs of the same or related traits, which were repeatedly detected in the different environments and/or populations . QTLs repeatedly detected in the same population from different environments (year-location combinations) were first integrated into identified QTLs, and then QTLs repeatedly detected in the different populations were integrated into consensus-QTLs. Each identified and consensus-QTLs were named according to a previous study .
SAM sampling and RNA isolation
To gain credible data, the SAMs were sampled from 10 representative individuals growing at the formation stage of floral bud primordia and then equally mixed for Zhongshuang11 and No.73290, respectively. Total RNA was isolated using an RNeasy Plant Mini Kit (Cat. 74124, Qiangen, Mississauga, ON) according to the manufacturer’s protocol. Then trace amount of DNA was treated with DNase I (Cat. 18068–015, Invitrogen, USA). In addition, the RNA quantity and quality were detected using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, USA), and the integrity of RNA was detected by electrophoresis with a 1% agarose gel in 1 ╳ TBE buffer at 70 V for 45 min.
cDNA libraries construction and Illumina sequencing
Sequencing libraries were prepared using NEB Next Ultra™ directional RNA Library Prep Kit for Illumina (San Diego, CA, USA) according to the manufacturer’s recommendations. Briefly, purification of mRNA was achieved from total RNA by using ploy-T oligo-attached magnetic beads. Then, adding the Illumina proprietary fragmentation buffer to the RNA samples, fragmentation of mRNA was completed under elevated temperature. Taking those fragments as templates, the first-strand cDNA was synthesized via using random hexamers and SuperScript II. The second-strand cDNA was synthesized by adding buffer, dNTPs, DNA polymerase I and RNase H to the reaction system. Subsequently, purification of the double-strand cDNA was performed using AMPure XP beads. The remaining overhangs were repaired via exonuclease/polymerase. Adenylation of the 3’ ends of cDNA fragments were conducted. Illumina PE adapter oligonucleotides were ligated to prepare for hybridization. 200-bp cDNA fragments were extracted using an AMPure XP system (Beckman Coulter, Beverly, CA, USA). DNA fragments were amplified by 12 cycles of PCR. After enrichment of DNA, the libraries were sequenced using Illumina Hiseq 2000 platform and raw data was generated. Adaptor, duplication and low quality sequences were removed to get clean data.
Differential gene expression analysis
First, all reads of each library were mapped to the reference genome using the SOAP program  with the default parameters. Second, the uniquely mapped reads were selected for quantifying the abundance. Third, the expression level was normalized using the values of RPKM (reads per kilobase per million reads). According to Bioconductor  package, DEseq (http://www.bioconductor.org/packages/release/bioc/html/DESeq.html) was used to measure gene differential expression between Zhongshuang11 and No.73290. The absolute value of log2 (Ratio) ≥ 1 (under the criterion of P ≤ 0.01 and false discovery rate (FDR) ≤ 0.001) was used as threshold to assess the significance of gene expression difference.
Functional annotation and category
Total DEGs were annotated based on BLSAT search in six public databases, including the protein family (Pfam) database (http://pfam.xfam.org/), InterPro database (http://www.ebi.ac.uk/interpro/entry/IPR001461), Swiss-Prot protein database (http://web.expasy.org/docs/swiss-prot_guideline.html), TrEMBLE database (http://web.expasy.org/docs/swiss-prot_guideline.html), the gene ontology (GO) database (http://geneontology.org/) and the Kyoto Encyclopedia of Genes and Genomes pathway (KEGG) database (http://www.genome.jp/kegg/pathway.html).
DEGs were functionally annotated using BLASTN with an E value cutoff of 1.0E-5 in TAIR 10 database. Then, The GO classification was performed using Classification SuperViewer Tool (http://bar.utoronto.ca/ntools/cgi-bin/ntools_classification_superviewer.cgi).
Pathway analysis was performed by searching against the Kyoto Encyclopedia of Gene and Genome (KEGG) pathway database with E value threshold of 1.0E-5.
Validation of RNA sequencing data by qRT-PCR
Twelve genes were selected randomly to validate RNA sequencing data by quantitative real-time PCR (qRT-PCR). The primer pairs were designed using Primer 5.0. Details of primer pairs were listed in Additional file 7: Table S7. Total RNA extracted for transcriptome sequencing was used for conducting qRT-PCR. qRT-PCR was performed using SYBR® Select Master Mix (2X) according to manufacturers’ recommendation. UBC21 gene was used as an internal control to normalize transcript levels . Real-time assay for each gene was performed with three independent biological replicates under identical conditions. Gene expression was calculated according to the previous study .
Differentially expressed gene
False discovery rate
Kyoto encyclopedia of genes and genomes
Quantitative trait loci
Recombinant inbred line
Reads per kilobase per million reads
Shoot apical meristem
Single nucleotide polymorphism
Diepenbrock W. Yield analysis of winter oilseed rape (Brassica napus L.): a review. Field Crop Res. 2000;67(1):35–49.
Li Y, Shen J, Wang T, Chen Q, Zhang X, Fu T, Meng J, Tu J, Ma C. QTL analysis of yield-related traits and their association with functional markers in Brassica napus L. Aust J Agric Res. 2007;58(8):759–66.
Shi J, Li R, Qiu D, Jiang C, Long Y, Morgan C, Bancroft I, Zhao J, Meng J. Unraveling the complex trait of crop yield with quantitative trait loci mapping in Brassica napus. Genetics. 2009;182(3):851–61.
Shi J, Zhan J, Yang Y, Ye J, Huang S, Li R, Wang X, Liu G, Wang H. Linkage and regional association analysis reveal two new tightly-linked major-QTLs for pod number and seed number per pod in rapeseed (Brassica napus L.). Sci Rep. 2015;5:14481.
Qu B, Zhang W, Chen XH, Li N, Cui N, Li TL. Research progress of flower bud differentiation mechanism of plant. Chin Agric Sci Bull. 2010;24:109–14.
Vaddepalli P, Scholz S, Schneitz K. Pattern formation during early floral development. Curr Opin Genet Develop. 2015;32:16–23.
Dodsworth S. A diverse and intricate signalling network regulates stem cell fate in the shoot apical meristem. Dev Biol. 2009;336(1):1–9.
Gaillochet C, Daum G, Lohmann JU. O cell, where art thou? The mechanisms of shoot meristem patterning. Curr Opin Plant Biol. 2015;23:91–7.
Holt AL, van Haperen JM, Groot EP, Laux T. Signaling in shoot and flower meristems of Arabidopsis thaliana. Curr Opin Plant Biol. 2014;17:96–102.
Zadnikova P, Simon R. How boundaries control plant development. Curr Opin Plant Biol. 2014;17:116–25.
Yruela I. Plant development regulation: Overview and perspectives. J Plant Physiol. 2015;182:62–78.
Suzaki T, Sato M, Ashikari M, Miyoshi M, Nagato Y, Hirano HY. The gene FLORAL ORGAN NUMBER1 regulates floral meristem size in rice and encodes a leucine-rich repeat receptor kinase orthologous to Arabidopsis CLAVATA1. Development. 2004;131(22):5649–57.
Suzaki T, Toriba T, Fujimoto M, Tsutsumi N, Kitano H, Hirano HY. Conservation and diversification of meristem maintenance mechanism in Oryza sativa: Function of the FLORAL ORGAN NUMBER2 gene. Plant Cell Physiol. 2006;47(12):1591–602.
Jiang L, Zhang W, Xia Z, Jiang G, Qian Q, Li A, Cheng Z, Zhu L, Mao L, Zhai W. A paracentric inversion suppresses genetic recombination at the FON3 locus with breakpoints corresponding to sequence gaps on rice chromosome 11L. Mol Genet Genomics. 2007;277(3):263–72.
Chu H, Qian Q, Liang W, Yin C, Tan H, Yao X, Yuan Z, Yang J, Huang H, Luo D, et al. The floral organ number4 gene encoding a putative ortholog of Arabidopsis CLAVATA3 regulates apical meristem size in rice. Plant Physiol. 2006;142(3):1039–52.
Xu J, Song X, Cheng Y, Zou X, Zeng L, Qiao X, Lu G, Fu G, Qu Z, Zhang X. Identification of QTLs for branch number in oilseed rape (Brassica napus L.). J Genet Genomics. 2014;41(10):557–9.
Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet. 2011;12(2):87–98.
Manuel G, Grabherr MG, Mitchell G, Cole T. Computational methods for transcriptome annotation and quantification using RNA-seq. Nat Methods. 2011;8(6):469–77.
Xu HM, Kong XD, Chen F, Huang JX, Lou XY, Zhao JY. Transcriptome analysis of Brassica napus pod using RNA-Seq and identification of lipid-related candidate genes. BMC Genomics. 2015;16(1):858.
Chalhoub B, Denoeud F, Liu S, Parkin IA, Tang H, Wang X, Chiquet J, Belcram H, Tong C, Samans B, et al. Plant genetics. Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome. Science. 2014;345(6199):950–3.
Li R, Li Y, Kristiansen K, Wang J. SOAP: short oligonucleotide alignment program. Bioinformatics. 2008;24(5):713–4.
Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. 2012;131(4):281–5.
Jiang Y, Deyholos MK. Transcriptome analysis of secondary-wall-enriched seed coat tissues of canola (Brassica napus L.). Plant Cell Rep. 2010;29(4):327–42.
Herrmann J, Ravilious GE, McKinney SE, Westfall CS, Lee SG, Baraniecka P, Giovannetti M, Kopriva S, Krishnan HB, Jez JM. Structure and mechanism of soybean ATP sulfurylase and the committed step in plant sulfur assimilation. J Biol Chem. 2014;289(15):10919–29.
Hao JH, Qi HY, Yan N, Wang HX. Advances in researches on flower bud differentiation of horticultural crops. Agric Sci Technol Equip. 2008;01:7–9.
Galston AW. Polyamines as Modulators of Plant Development. Bioscience. 1983;33(6):382–8.
Sun WQ, Chu MY. Study on content variations of endogenous amino acids at the critical period of floral initiation in MEI-FLOWER. Acta Agriculiurae Shanghai. 1990;2:43–7.
Sassi M, Ali O, Boudon F, Cloarec G, Abad U, Cellier C, Chen X, Gilles B, Milani P, Friml J, et al. An auxin-mediated shift toward growth isotropy promotes organ formation at the shoot meristem in Arabidopsis. Curr Biol. 2014;24(19):2335–42.
Wagner D, Sablowski RW, Meyerowitz EM. Transcriptional activation of APETALA1 by LEAFY. Science. 1999;285(5427):582–4.
Weigel D, Alvarez J, Smyth DR, Yanofsky MF, Meyerowitz EM. LEAFY controls floral meristem identity in Arabidopsis. Cell. 1992;69(5):843–59.
Baumann K, Venail J, Berbel A, Domenech MJ, Money T, Conti L, Hanzawa Y, Madueno F, Bradley D. Changing the spatial pattern of TFL1 expression reveals its key role in the shoot meristem in controlling Arabidopsis flowering architecture. J Exp Bot. 2015;66(15):4769–80.
Takacs EM, Li J, Du C, Ponnala L, Janick-Buckner D, Yu J, Muehlbauer GJ, Schnable PS, Timmermans MCP, Sun Q, et al. Ontogeny of the Maize Shoot Apical Meristem. Plant Cell. 2012;24(8):3219–34.
Wong CE, Singh MB, Bhalla PL. The dynamics of soybean leaf and shoot apical meristem transcriptome undergoing floral initiation process. PLoS One. 2013;8(6):e65319.
Murray JA, Jones A, Godin C, Traas J. Systems analysis of shoot apical meristem growth and development: integrating hormonal and mechanical signaling. Plant Cell. 2012;24(10):3907–19.
Shani E, Yanai O, Ori N. The role of hormones in shoot apical meristem function. Curr Opin Plant Biol. 2006;9(5):484–9.
Gendron JM, Jiang-Shu L, Min F, Ming-Yi B, Stephan W, Springer PS, Kathryn B, Zhi-Yong W. Brassinosteroids regulate organ boundary formation in the shoot apical meristem of Arabidopsis. Proc Natl Acad Sci U S A. 2012;109(51):21152–7.
Han Y, Gao S, Muegge K, Zhang W, Zhou B. Advanced Applications of RNA Sequencing and Challenges. Bioinformatics Biol Insights. 2015;9 Suppl 1:29–46.
Wang Z, Gerstein MM. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.
Fernandez-Lozano A, Yuste-Lisbona FJ, Perez-Martin F, Pineda B, Moreno V, Lozano R, Angosto T. Mutation at the tomato excessive number of floral organs (ENO) locus impairs floral meristem development, thus promoting an increased number of floral organs and fruit size. Plant Sci. 2015;232:41–8.
Sun B, Ito T. Regulation of floral stem cell termination in Arabidopsis. Front Plant Sci. 2015;6:17.
Han Y, Jiao Y. APETALA1 establishes determinate floral meristem through regulating cytokinins homeostasis in Arabidopsis. Plant Signal Beh. 2015;10(11):e989039.
Bowman JL, Alvarez J, Weigel D, Meyerowitz EM, Smyth DR. Control of flower development in Arabidopsis thaliana by APETALA1 and interacting genes. Development. 1991;119(3):721–43.
Irish VF, Sussex IM. Function of the apetala-1 gene during Arabidopsis floral development. Plant Cell. 1990;2(8):741–53.
Andres F, Romera-Branchat M, Martinez-Gallegos R, Patel V, Schneeberger K, Jang S, Altmuller J, Nurnberg P, Coupland G. Floral Induction in Arabidopsis by FLOWERING LOCUS T Requires Direct Repression of BLADE-ON-PETIOLE Genes by the Homeodomain Protein PENNYWISE. Plant Physiol. 2015;169(3):2187–99.
Khan M, Ragni L, Tabb P, Salasini BC, Chatfield S, Datla R, Lock J, Kuai X, Despres C, Proveniers M, et al. Repression of Lateral Organ Boundary Genes by PENNYWISE and POUND-FOOLISH Is Essential for Meristem Maintenance and Flowering in Arabidopsis. Plant Physiol. 2015;169(3):2166–86.
Ooka H, Satoh K, Doi K, Nagata T, Otomo Y, Murakami K, Matsubara K, Osato N, Kawai J, Carninci P, et al. Comprehensive analysis of NAC family genes in Oryza sativa and Arabidopsis thaliana. DNA Res. 2003;10(6):239–47.
Tromas A, Paque S, Stierle V, Quettier AL, Muller P, Lechner E, Genschik P, Perrot-Rechenmann C. Auxin-binding protein 1 is a negative regulator of the SCF(TIR1/AFB) pathway. Nat Commun. 2013;4:2496.
Dornelas MC, Van Lammeren AA, Kreis M. Arabidopsis thaliana SHAGGY-related protein kinases (AtSK11 and 12) function in perianth and gynoecium development. Plant J. 2000;21(5):419–29.
Shuai B, Reynaga-Pena CG, Springer PS. The lateral organ boundaries gene defines a novel, plant-specific gene family. Plant Physiol. 2002;129(2):747–61.
Nakazawa M, Ichikawa T, Ishikawa A, Kobayashi H, Tsuhara Y, Kawashima M, Suzuki K, Muto S, Matsui M. Activation tagging, a novel tool to dissect the functions of a gene family. Plant J. 2003;34(5):741–50.
Ho WW, Weigel D. Structural features determining flower-promoting activity of Arabidopsis FLOWERING LOCUS T. Plant Cell. 2014;26(2):552–64.
Ping J, Liu Y, Sun L, Zhao M, Li Y, She M, Sui Y, Lin F, Liu X, Tang Z, et al. Dt2Is a Gain-of-Function MADS-Domain Factor Gene That Specifies Semideterminacy in Soybean. Plant Cell. 2014;26(7):2831–42.
Randoux M, Daviere JM, Jeauffre J, Thouroude T, Pierre S, Toualbia Y, Perrotte J, Reynoird JP, Jammes MJ, Hibrand-Saint Oyant L, et al. RoKSN, a floral repressor, forms protein complexes with RoFD and RoFT to regulate vegetative and reproductive development in rose. New Phytol. 2014;202(1):161–73.
Wang Y, Pijut PM. Isolation and characterization of a TERMINAL FLOWER 1 homolog from Prunus serotina Ehrh. Tree Physiol. 2013;33(8):855–65.
Wang X, Xie B, Zhu M, Zhang Z, Hong Z. Nucleostemin-like 1 is required for embryogenesis and leaf development in Arabidopsis. Plant Mol Biol. 2011;78(1–2):31–44.
Wang X, Gingrich DK, Deng Y, Hong Z. A nucleostemin-like GTPase required for normal apical and floral meristem development in Arabidopsis. Mol Biol Cell. 2012;23(8):1446–56.
Veit B, Briggs SP, Schmidt RJ, Yanofsky MF, Hake S. Regulation of leaf initiation by the terminal ear 1 gene of maize. Nature. 1998;393(6681):166–8.
Doyle JJ. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem Bull. 1987;19:11–5.
Goffinet B, Gerber S. Quantitative trait loci: a meta-analysis. Genetics. 2000;155(1):463–73.
Li N, Shi J, Wang X, Liu G, Wang H. A combined linkage and regional association mapping validation and fine mapping of two major pleiotropic QTLs for seed weight and silique length in rapeseed (Brassica napus L.). BMC Plant Biol. 2014;14:114.
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5(10):R80.
Chen X, Truksa M, Shah S, Weselake RJ. A survey of quantitative real-time polymerase chain reaction internal reference genes for expression studies in Brassica napus. Anal Biochem. 2010;405(1):138–40.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) Method. Methods. 2001;25(4):402–8.
Part of the filed work was carried out at the farm of Henan Academy of Agricultural Sciences. The authors thank Professor Shufen Zhang, Dr. Junping He and Mrs Yujiao Liu for their kind help in field experiment and management.
This research was supported by the National Basic Research and Development Program (2015CB150202), the Rapeseed Industry Technology System (CARS-13), the Agricultural Science and Technology Innovation Project (CAAS-ASTIP-2013-OCRI), the Core Research Budget of the Non-profit Governmental Research Institution (1610172014001), and the Hubei Agricultural Science and Technology Innovation Center of China.
Availability of data and materials
The datasets supporting the conclusions of this article are available in the Sequence Read Archive at the National Center for Biotechnology Information (NCBI) at the accession number SRX1715587. The datasets supporting the conclusions of this article are included within the article and its additional files.
JS and HW designed the experiments; JS performed the construction, genotyping and phenotyping of the RIL population; YY collected samples and prepared RNA. JY performed QTL scanning and expression analysis. XW and GL performed the field experiments; JY and JS conducted the data analysis and wrote the manuscript. All the authors read and approved the final manuscript.
The authors declare no competing financial interests.
Consent for publication
Ethics approval and consent to participate
ANOVA analysis of pod number in four experiments. PNm/PNb/PNw was the abbreviation of pod number from the main inflorescence, branch inflorescence and whole plant (XLS 28 kb)
A. List of all detected QTLs at both significant and suggestive level for pod number. B. List of identified QTLs after deleting non-reproducible suggestive QTLs for pod number. C. Details of identified QTLs for pod number in the BnaZNRIL population in rapeseed. (XLSX 21 kb)
List of QTLs identified in the our and other studies for pod number in rapeseed. (XLSX 21 kb)
The differentially expressed genes (DEGs) between Zhongshuang11 and No.73290 samples using false discovery rate (FDR) < 0.05 and fold change > 2 as the significant cutoff. RPKM: reads per kilobase per million reads. (XLSX 1029 kb)
Annotation for DEGs (E value < 1E-5). The most matched one was used to annotate the rapeseed genes. (XLSX 218 kb)
190 KEGG pathways represented by DEGs between Zhongshuang11 and No.73290. (XLSX 128 kb)
Gene-specific primers for quantitative real-time PCR. UBC21 is the internal reference gene. F and R denote forward and reverse primers. (XLSX 10 kb)
About this article
Cite this article
Ye, J., Yang, Y., Chen, B. et al. An integrated analysis of QTL mapping and RNA sequencing provides further insights and promising candidates for pod number variation in rapeseed (Brassica napus L.). BMC Genomics 18, 71 (2017). https://doi.org/10.1186/s12864-016-3402-y
- Brassica napus
- Pod number
- QTL mapping
- RNA sequencing
- Candidate gene