Exploration into Natural Variation for Genes Associated with Determinate and Capitulum-like Inflorescence in Brassica napus

1. Introduction

Rapeseed (Brassica napus L., AACC, 2n=38) is the world’ third most important source of vegetable oil after palm and soybean [1]. In China, rapeseed accounts for about 85% acreage of oilseed Brassica, and provided the second most important vegetable oil after soybean [2]. The great importance makes it an ideal model species for theory and application research.

The precise development of inflorescences and flowers is crucial for reproductive success in flowering plants. B. napus is a simple structure typical of the Brassicaceae, and develops an indeterminate, racemetype inflorescence comprising individual lateral flowers arising immediately and sequentially from an apical inflorescence meristem (IM) [3,4]. The lateral flowers in inflorescence developed according to a well-defined program of events that gives rise to a stereotypical floral structure comprising a fixed sequence of concentric whorls with fixed numbers of floral organs (four sepals, four petals, six stamens and a central pistil) [4]. While once the flowers arise on the top of IM, naming terminal flower, the growth of inflorescence will be ceased, which is defined as determinate inflorescence. The situation of terminal flower had been found in many species, such as in Arabidopsis thaliana [5], Glycine max [6], Nicotiana tabacum [7], Brassica juncea [8], Sesamum indicum [9], Brassica napus [10].

The characteristics of inflorescences and flowers were involved in one complex genetic network, among which LEAFY (LFY), APETALA1 (AP1) and TERMINAL FLOWER1 (TFL1) played one key role. LFY is one key transcriptional regulator in the network establishing flower initiation in floral meristem (FM), and is activated by genes AGL24, SOC1, MP and ANT/AIL6 [11,12,13,14,15]. For the establishment of floral meristem identity, AP1 functions a key role and is activated by LFY, CAL and LMI2 [16]. The TFL1 gene had been drawn a widely interests for its important role in shifting indeterminate inflorescence to determinate inflorescence due to the appearance of terminal flower in the top of IM. The terminal flower mutant was firstly isolated from the recessive mutations by screening a M₂ population derived from EMS mutagenized seeds of Arabidopsis thaliana ecotype Columbia [5]. Then the candidate recessive gene of tfl1 for regulating the terminal flower phenotype in Arabidopsis thaliana was mapped on the top arm of chromosome 5 [17] and cloned [18]. After that, tfl1′ homolog in Nicotiana tabacum [7] and Glycine max [6] were cloned. In Sesamum indicum, the Sidt1 gene homologous to Arabidopsis TFL1 was mapped on LG09 in a genome region of 41 kb by one ultra-dense SNP genetic map [19]. The Sdt1 homologous to Arabidopsis TFL1 gene associated with determinate feature in Brassica juncea was mapped to the linkage group B05, which was flanked by SSR markers SJ6842 and Ni4-A10 at distances of 15.9 cM and 14.0 cM, respectively [8]. One recent paper in B. napus published the discovery of one microspore culture-origin determinate mutation with terminal flower [20]. The regulator Bnsdt1, homologous to Arabidopsis TFL1, was fine mapped on one region of approximately 220 kb, between 16,627 and 16,847 kb on A10 using BC1 and BC3 populations [20]. The cooperative function of these genes could result in differential inflorescence architectures. “A unifying inflorescence model” postulated that TFL1 could increase vegetativeness (veg) and LFY could reduce veg in meristems, leading to different architectures in Arabidopsis thaliana [3]. In another model, increased TFL1 expression could lead to larger inflorescences with more and longer branches, whereas increased AP1 expression could lead to smaller inflorescences with fewer branches and flowers [21].

Except for the traditional techniques for genetic research, such as QTL mapping, many novel techniques based on sequencing were numerously emerging in recent years. RNA-Seq is a recently developed approach to transcriptome profiling using deep-sequencing technologies [22]. The obtained transcriptome included all coding mRNA and noncoding RNA sequences in the cells of one specific development stage or physiological condition. Plant transcriptome analysis is fast and considerable for providing information on highly expressed genes, differentially expressed genes, new genes for function analysis and gene screening related to studied traits [23,24,25,26]. Besides, BSA-Seq strategy was one more efficient QTL mapping method comparing with the traditional QTL mapping method. Combing the traditional BSA method and next-generation sequencing (NGS) technique, BSA-Seq technique has been widely used in QTL mapping for the precise identification of target genes [27,28,29], such as in rapeseed [30], rice [29], maize [31] and so on. Furthermore, the combined analysis of BSA-seq and RNA-seq data enabled the identification of candidate genes [32,33,34].

In the present study, we reported one natural mutant (MT) with determinate and capitulum-like inflorescence in B. napus. Then phenotypic analysis, genetic analysis, RNA-Seq and BSA-Seq methods were used to identify the candidate genes for regulating the determinate inflorescence. The results would be helpful in studying the molecular mechanism of the mutation phenotype and improving the inflorescence architecture in B. napus.

2. Results

2.1. Phenotypic characterization of the mutant (MT) plants

One mutant (MT) plant (Figure 1B/b and Figure 1D/d) was found in the self-polinated progenies of the wide-type (WT) B. napus line of GD605-2 (Figure 1A/a and Figure 1C/c) in the field. The phenotype of MT plants and the WT plants was kept the same until the plants started to bolt. Then at the apex of each inflorescence axis, the indeterminate inflorescence was terminated with a terminal flower in MT plants (Figure 1B/b and Figure 1D/d). The terminal flowers contained one pistil and six stamens, while no petals and seeds would be formed. Each inflorescence axis with buds/flowers at the apex closer to the “terminal flower” was kept to be clustered together and not prolongated, which finally developed into one capitulum-like head in MT plants (Figure 1A/a and Figure 1D/d). Each capitulum-like head had about 10.0±0.9 normal flowers/siliques in the top of each inflorescence axis. The MT plants had a plant height of about 158.8±4.3 cm, significantly (p=0.001) lower than the WT plants of 180.9±3.4 cm. Besides, no significant difference (p=0.903) in the main inflorescence was detected for the number of siliques between WT plants (58.3±8.2) and MT plants (58.0±9.9).

The development of inflorescence apex from three to seven-leaf stages of MT and WT plants was observed by paraffin sectioning technique (Figure 2). The inflorescence apex from three to four-leaf stages developed normally in all the detected MT and WT plants (Figure 2A/a and Figure 2B/b). Starting from five-leaf stage, the inflorescence apex began to split and developed into one terminal flower (Figures D-F/d-f).

2.2. Genetic analysis of the MT phenotype

Genetic determinism behind the MT phenotype was analyzed through the phenotypic characterization of F₁ progenies from crosses between the MT and WT plants. All the F₁ plants from 4 reciprocal crosses between MT and Ningyou7 (29 and 33 plants, respectively), Bakow (29 and 33 plants, respectively), 97009 (26 and 30 plants, respectively), 97081 (28 and 15 plants, respectively) displayed similar phenotype of indeterminate inflorescences to the WT parent plants. Since all the F₁ plants displayed the normal phenotype of indeterminate inflorescences, the allele(s) responsible for the MT phenotype showed complete recessiveness over the allele(s) responsible for the normal phenotype of indeterminate inflorescences.

In the F₂ population from MT×WT, 553 plants with indeterminate inflorescences and 131 plants with determinate inflorescences in 2018, fitting with a segregation ratio of 13:3 (χ²=0.073, p= 0.788). In 2020, the F₂ population from MT×WT was classified into two groups: plants with indeterminate inflorescences (1,077 plants) and plants with determinate inflorescences (236 plants), fitting with a segregation ratio of 13:3 (χ²=0.519, p= 0.471). Furthermore, another F₂ population from WT×MT in 2020 was also showed the same segregation ratio of 13:3 (843 plants with indeterminate inflorescences versus 172 plants with determinate inflorescences (χ²=2.169, p= 0.141). Therefore, the segregation patterns of inflorescence characterization in the F₂ populations supported a digenic inheritance model in these crosss.

2.3. QTL mapping by BSA-Seq technique

BSA-Seq technique was used to rapidly map the QTL accounting for the determinate inflorescences of MT plants. The “determinate” DNA pool, “indeterminate” DNA pool and the two parental DNA pools were subjected to Illumina sequencing with an average sequencing depths of 16.25-fold. Then the clean reads were aligned to the reference genome by using BWA software, and the average proportion of mapped reads to clean reads was 97.08%. After removing low quality reads, a total of 94.04 Gbp clean data were obtained with an average Q30 rate of 80% and average of onefold coverage ratio of 92.16%. All the obtained clean data were used to develop SNP and InDel. A total of 1,500,703 SNPs and 945,276 SNPs were identified between two parents and 945,276 ones between the two DNA pools, respectively. Furthermore, a total of 400,329 InDel between two parents and 280,916 InDel between two DNA pools were identified. The ΔSNP-index and ΔInDel-index were used for QTL mapping of the candidate intervals for determinate inflorescence. As a result, two overlabbed regions were identified on C02, wihch had a size of 4.18 Mb (14.27–18.45 Mb) byΔSNP-index and 4.18 Mb (14.23–18.41 Mb) byΔInDel-index (Figure 3), respectively. On C06, two another overlabbed regions were identified on C06, wihch had a size of 0.7 Mb (32.98–33.68 Mb) byΔSNP-index and 0.72 Mb (32.97–33.69 Mb) byΔInDel-index, respectively. Finally, these two overlapped regions on C02 (14.27–18.41 Mb, 320 candidate genes) and C06 (32.98–33.68 Mb, 117 candidate genes) were chosen as the candidate regions for genetic control of determinate inflorescences in MT plants.

2.4. Gene expression profile analysis

A total of six cDNA libraries were constructed and sequenced, an average of 48,705,785 and 46,319,327 raw reads was generated from the WT and MT libraries, respectively (Table 1 and Table S1). After removing low quality reads, adapter polluted reads and higher N content (>5%) reads, an average of 47,184,448 (WT) and 44,944,151 (MT) clean reads were obtained. After blasting the reference genome of B. napus and filtering the genes that only contained one exon or encoded short peptide chains (<50 amino acid residues), a total of 86,026 genes were revealed through blasting the reference genome using DESeq2 (v1.6.3). It was found that 6,309 genes were differentially expressed, among which 4,103 ones were down-regulated and 2,209 ones were up-regulated in MT compared with WT (Table S2). Based on the FPKM values, the up-regulated and down-regulated genes between WT and MT were also revealed by a hierarchical clustering analysis in Figure S1. To verify the expression of DEGs detected by RNA-Seq, 16 candidate genes (DEGs) regulating the MT phenotype were randomly chosen for validation by qRT-PCR. The data obtained by qRT-PCR was consistent with the RNA-Seq results (Figure S2 and Table S3), suggesting the reliability of the transcriptome database.

To functionally annotate the B. napus transcriptome, the 6,309 DEGs were blasted in search of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Finally, 4,548 DEGs were successfully annotated, among which 4,522 (71.68%) in GO and 2,127 (33.71%) in KEGG (Table S4). For the GO classification analysis of DEGs, 4,522 genes were assigned into three main GO functional categories (cellular component, biological process and molecular function) and then divided into 55 sub-categories (Figure S3 and Table S5). For the KEGG analysis of gene functions, all the annotated 6,309 DEGs were assigned into 128 pathways (Table S6) based on KEGG database. Several DEGs were assigned into more than one sub-category.

2.5. Candidate genes for regulating the determinate inflorescence

Following the annotation of GO and KEGG, the potential DEGs for regulating the flower and inflorescence mutation were further blasted with the database of Arabidopsis thaliana (https://www.arabidopsis.org/index.jsp) and NCBI (https://www.ncbi.nlm.nih.gov/) for detailed information. A total of 133 candidate genes for regulating the flower development (75 genes, 58.6%), shoot meristem development (29 genes, 22.7%) and inflorescence meristem development were identified (13 genes, 9.8%) (Figure. S4 and Table S7). Several genes were involved in more than one sub-category.

To integrate the results of BSA-Seq and RNA-sequencing, we perform an alignment analysis between the 133 candidate genes potentially related to determinate inflorescence of the MT plants and the reference genome of B. napus (https://www.genoscope.cns.fr/brassicanapus/) by the BLAST-like alignment tool [35]. Eight DEGs corresponding to seven genes (Table 2) located on C02 between 14.27 Mb and 18.41 Mb were detected, among which the widely recognized genes for regulating determinate inflorescence of BnaC02g02900D encoding TERMINAL FLOWER 1 (TFL1) (ChrC02: 1,320,657–1,321,719) was included [20,36]. In addition, three DEGs on C06 (Table 2) between 32.98 Mb and 33.68 Mb were detected, among which only one DEG of BnaC06g25500D encoding the gene of APETALA1 (AP1) (C06: 27,150,336–27,153,999) was related to the inflorescence development. Then the two candidate genes were renamed as BncTFL1 and BncAP1. The expression level of BncTFL1 and BncAP1 showed significantly increased from 5-leaf stage and started to decrease from 7-leaf stage in WT plants (Figure 4). In addition, the expression level of BncTFL1 and BncAP1 was significantly decreased from 5-leaf stage to 7-leaf stage in MT plants comparing with WT plants (Figure 4).

3. Discussion

In the present study, one natural mutation with determinate inflorescences characterizing with terminal flower and capitulum-like inflorescences in Brassiceae were reported. The mutation was characterized by paraffin sectioning, and the inflorescence apex began to split and developed into one terminal flower, starting from five-leaf stage. The segregation patterns of inflorescences in the F₂ populations supported a digenic inheritance model, which was further approved by BSA-Seq technique. TERMINAL FLOWER 1 (TFL1) had been proved to be responsible for controlling the determinate inflorescences in many species [5,6,7,8,9,20]. While BncTFL1 mapped on C02 in the present study was one new allele of TFL1, differently from these genes on B05 of B. juncea [8] and on A10 of B. napus [10] in Brassiceae. In addition, AP1 had been identified as the genes controlling the development of inflorescence meristem and floral meristem [21], BncAP1 on C06 was first reported as the candidate gene of controlling determinate inflorescences. Quite different from the reported phenotype of determinate inflorescences, the capitulum-like inflorescences was also found in GD605-2. The formation of capitulum-like inflorescences had once been well exploited by two models of FL1 & LFY [3] and TFL1 & AP1 [21], in which the specific combination of gene expression levels result in the emergence of different capitulum-like inflorescences.

The special mutation with determinate and capitulum-like inflorescence could provide unique materials for the basic research of the development of inflorescence in B. napus. To broaden the utilization of the MT phenotype, we had transferred the mutation phenotype into germplasm resources with colored flowers, purple stem and leaves, muti-head inflorescence. These new materials could provide useful resources for rapeseed tourism. In the future, we would focus on the gene cloning and functional analysis of these two candidate genes of BncTFL1 and BncAP1.

4. Materials and methods

4.1. Plant materials and field experiment

In March 2014, one mutant (MT) plant with determinate and capitulum-like inflorescence was found. The MT plant was naturally mutated and found in the self-polinated progenies of one breeding line GD605-2 (F₇) (wild-type, WT). In October 2017, twenty plants for each of the MT (F₄) and WT (F₇) lines were planted in the field of Guizhou University, Guiyang, China. In March 2018, the inflorescences from three randomly chosen plants for each of MT and WT lines were sampled, immediately frozen in liquid nitrogen, and stored at -80⁰C for RNA extraction.

4.2. Investigation of the agronomic traits of MT and WT plants

In April 2019, twenty plants for each of MT and WT plants were used for investigation of number of flowers/siliques in the capitulum-like head, plant height, number of siliques in the main inflorescence.

4.3. Paraffin sectioning

The paraffin sectioning is mainly following the work of Zhou, et al. [37] with a minor modification. The shoot tips of MT and WT lines from three to seven-leaf stages are dissected as fast as possible after sampling. The dissected tissue samples were fixed for >24 h in FAA fixative in chemical hood. The fixed tissues were then continuously dehydrated and cleared, respectively, using ethanol and xylene in volume ratios of 3:1, 1:1, and 1:3. After the tissue samples were in melted paraffin, they were removed from the mold and embedded in 60 °C paraffin wax solution which was allowed to solidify at room temperature. Paraffin sections (10 μm) were obtained using an automatic microtome. The sections were stained with aniline blue, and observed with a microscope. Images of the anthers at different stages were captured with a camera.

4.4. BSA-Seq analysis

The BSA-Seq method was used for mapping the genes of regulating the phenotype of determinate inflorescence. The “determinate” bulk was made by mixing equal amount of DNA from 10 F₂ plants with determinate inflorescence, while the “indeterminate” bulk was formed from 10 F₂ plants with indeterminate inflorescence. The two bulks and the DNA samples of the parental lines of WT plants and MT plants were sequenced on an Illumina HiSeq™2000 platform (Beijing Biomarker Biotechnology Co., Beijing, China). The low-quality reads containing adaptors were filtered. The reads with more than 10% of missing bases and more than 50% of bases with Q-score lower than 10 were filtered, and the clean reads thus obtained were mapped to the Brassica napus reference genome (Genoscope v4.1, http://www.genoscope.cns.fr/brassicanapus/data/) using BWA software [29,38]. SNP-calling and SNP annotation were done using the professional softwares of ANNOVAR and SAM [38,39], MarkDuplicatePicards (http://sourceforge.net/projects/picard/), GATK [40], SnpEff [41]. The differences in allele frequency between bulked pools was used to calculate the SNP-index for identifying the candidate regions of the genome associated with determinate inflorescences [28,29,42]. In detial, Δ(SNP index) was calculated by sliding window analysis among the genome within 1 Mb width windows and 1 kb at each step [28]. Totally, all those analysis was performed with the related tools on the online open platform of BMKCloud (http://www.biocloud.com/).

4.5. RNA extraction, preparation, sequencing and data analysis

Total RNA (2 µg) was extracted from the shoots (0.2–0.5g) of three independent plants for each of WT and MT lines in 5-leaf stage using the TRIzol kit (Invitrogen, Carlsbad, CA), according to the manufacturer’s instructions. The RNA purity was checked using the Kaiao K5500^®Spectrophotometer (Kaiao, Beijing, China), and the RNA integrity and concentration were assessed using the RNA Nano 6000 Assay Kit for the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Then, the six RNA samples were sent to the ANOROAD GENOME company (http://www.genome.cn/) for the construction of cDNA libraries and Illumina deep sequencing according to the paper of Wang, et al. [43]. The raw RNA-sequencing data were filtered by a Perl script, following the steps of Wu, et al. [44] and our earlier study [34].

4.6. Identification and annotation of differentially expressed genes (DEGs)

DESeq2 v1.6.3 was used for differential gene expression analysis between MT and WT with three biological replicates under the theoretical basis obeys the hypothesis of negative binomial distribution for the value of count. The p-value was corrected by the BH method. Genes with q≤0.05 and |log2_ratio|≥1 were identified as differentially expressed genes (DEGs) [45]. The DEGs obtained were further annotated with Gene Ontology (GO, http://geneontology.org/) and analyzed by KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.kegg.jp/) [46,47]. The GO enrichment of DEGs was implemented by the hypergeometric test, in which the p-value is calculated and adjusted to produce the q-value, and the data background is the genes in the whole genome. GO terms with q<0.05 were considered to be significantly enriched. GO enrichment analysis was used to determine the biological functions of the DEGs. The KEGG enrichment of the DEGs was determined by the hypergeometric test, in which p-value was adjusted by multiple comparisons to produce the q-value. KEGG terms with q<0.05 were considered to be significantly enriched.

4.7. Quantitative real time-PCR (qRT-PCR) analysis

Quantitative real-time PCR (qRT-PCR) was used to verify the transcript levels of the RNA-Seq results. Total RNA was extracted using the TRIzol kit (Invitrogen), according to the manufacturer’s instructions. Then, the cDNA was synthesized by reverse transcription using PrimeScript RT reagent kits with gDNA Eraser (Takara, Dalian, China) according to the manufacturer’s instructions. Sixteen gene-specific primers for qRT-PCR were designed based on reference unigene sequences randomly chosen from the DEGs using Primer Premier 5.0. Real-time PCR was conducted using SsoAdvancedTM Universal SYBR Green Supermix (Hercules, CA) according to our earlier research [34]. The 2-ΔΔCt algorithm was used to calculate the relative level of gene expression. The β-actin gene was used as the internal control, and the WT samples served as the control. All qRT-PCR were performed with three biological replicates, and run on a Bio-Rad CFX96 Real Time System (Bio-Rad, Hercules, CA, USA).

5. Conclusions

In the present study, we reported one natural mutation with determinate and capitulum-like inflorescence in B. napus. The mutation started from 5-leaf stage in the shoots to the flowering stage in inflorescence with a phenotype of one top flower and capitulum-like inflorescence. Trough genetic analysis and BSA-Seq analysis, two QTL regions on C02 and C06 were detected. To better dig the candidate genes within the QTL regions, RNA-Seq analysis was done. Finally, one joint analysis combing BSA-Seq and RNA-Seq identified two candidate gene of BncTFL1 and BncAP1 for regulating the MT phenotype. Except for one useful material for basic research of inflorescence development, the mutation could also be used for rapeseed tourism as for its special inflorescence phenotype.