Investigation for a multi-silique trait in Brassica napus by alternative splicing analysis

Introduction

Rapeseed (Brassica napus L.), an allotetraploid with a complex genome (AACC, 2n = 38), is the second leading source of vegetable oil globally (Liu et al., 2015). The agronomic traits related to rapeseed yield include the pod (silique) number per plant, branch number, and seed weight (Liu et al., 2015; Zhang et al., 2006; Li et al., 2015). We previously reported that zws-ms, a multi-silique rapeseed line (Chai et al., 2019), produces three independent pistils and 9 to 10 stamens on the same receptacle in a flower, which consequently leads to the formation of three independent siliques on a carpopodium rather than the single siliques typically observed. Moreover, this trait was found to be affected by the environment, with temperature considered to be the factor most likely to switch on/off the formation of multi-silique.

Temperature is a major environmental factor that regulates various aspects of plant morphology, physiology, and biochemistry, affecting germination, growth, development, and flowering (Ren et al., 2019). Fertility in crops such as rapeseed (Yu et al., 2015) and rice (Oryza sativa) (Yu et al., 2017) is affected by temperature. In winter rapeseed lines, although a period of vernalization under low temperature is necessary to initiate flowering, cold stress inhibits growth and development, disturbs metabolism, and causes wilting or even death. Notably, cold stress also induces alternative splicing (AS) in plants (Palusa, Ali & Reddy, 2007; Iida, 2004).

AS is defined as the mechanism by which primary transcripts are processed into two or more mature isoforms, which enables a single gene to produce diverse protein products (Pan et al., 2008; Sablok et al., 2011). These proteins differ from each other not only in structure but also possibly in function, subcellular localization, and/or stability (Huang et al., 2019; Chauhan et al., 2019). AS is common in plants; for example, in Arabidopsis thaliana, more than 60% of intron-containing genes undergo AS (Syed et al., 2012). Many environmental factors regulate AS events in plants, including CO₂ concentration (Huang et al., 2019), light (Godoy et al., 2019), salt stress (Ding et al., 2014), and nutrient deficiencies (Nishida et al., 2017). AS not only provides an important source of transcriptomic and proteomic diversity and plasticity for use in natural selection (Labadorf et al., 2010), but it also plays specific roles in the response (Chauhan et al., 2019) or adaptation to environmental stresses (Filichkin et al., 2015). Guo et al. (2019a) and Guo et al. (2019b) identified four splicing variants of two BnCYCD3-1-LIKE genes in B. napus and found evidence that their AS may play an important role in the response to environmental stresses. Xia et al. (2017) discovered that the AS with intron retention of EARLY MATURITY8 (EAM8) led to early flowering in a barley (Hordeum vulgare) landrace; while in shepherd’s purse (Capsella bursa-pastoris), flowering time varied with changes in the splicing of a FLOWERING LOCUS C (FLC) homolog (Slotte et al., 2009). In addition, the heterologous expression of a vacuolar membrane Na⁺/H⁺ antiporter gene (SsNHX1) AS variant from seepweed (Suaeda salsa) enhanced the salt tolerance of Arabidopsis (Li et al., 2009).

As mentioned above, low temperatures switch off the multi-silique trait in zws-ms rapeseed. When zws-ms plants were planted in Xindu, Sichuan Province, China, the multi-silique trait was continuously stable for years; however, when they were grown in Ma’erkang, Sichuan Province, where the annual average temperature is consistently 7.6 °C lower, the multi-silique trait disappeared and all plants displayed normal siliques (Chai et al., 2019). We previously investigated the association of chromosomal regions with this trait, at the genomic and transcriptomic levels, selecting potential candidates from the differentially expressed genes (DEGs) between the multi- and single-silique plants. However, the involvement of post-transcriptional modifications and the mechanisms by which temperature regulates this multi-silique trait remain unclear. AS is often responsive to cold stress in plants (Iida, 2004; Palusa, Ali & Reddy, 2007) and is a mechanism by which plants perceive temperature fluctuations and modulate the activity of their transcription factors (Seo, Park & Park, 2013). In view of the above insights, we analyzed AS using transcriptome sequencing (RNA-seq) in this study. High-throughput RNA-seq technology is a widely used, highly efficient, and is an economical strategy for transcriptomic profiling (Tong et al., 2013; Wang, Gerstein & Snyder, 2009). It has become increasingly popular because of the following qualities (Mortazavi et al., 2008; Ozsolak & Milos, 2011); (Marioni et al., 2008; Tong et al., 2013; Sultan et al., 2008): (1) it can be used to detect and quantify the expression of genes, including those expressed at low levels; (2) it can facilitate the annotation of genes and lead to the discovery of novel genes or transcripts; (3) the results are highly reproducible between both technical and biological replicates; and (4) it can detect AS events.

We performed transcriptome sequencing (RNA-seq) on the flower buds of zws-ms and its near-isogenic line (NIL), zws-217, which produces normal single siliques. This facilitated the identification of the AS events in both lines and the analysis of the differentially expressed AS genes and those with line-specific AS events. Combining these data with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations, we identified likely candidate genes switching on or off the multi-silique trait by altering AS events or transcriptional levels in varied environments.. To the best of our knowledge, this is the first time that the regulation of flower/fruit morphology by AS has been investigated in rapeseed, and our results provide insights into this field more generally.

Materials & Methods

Plant materials and growth conditions

The rapeseed line zws-ms and its NIL, zws-217 (Chai et al., 2019), were kept in the Crop Research Institute, Sichuan Academy of Agricultural Sciences, China. Both zws-ms and zws-217 were homozygous for almost all genes, differing from each other only in the multi-silique trait of zws-ms (Fig. 1). The NILs zws-217 and zws-ms were both grown in an experimental field in the Xindu District of Chengdu in the Sichuan Basin, China, under normal environmental conditions. Additionally, the both lines were also grown in Ma’erkang, a mountainous area in western Sichuan, with a much lower annual average temperature. The annual average temperature in Xindu and Ma’erkang is 16.2 °C and 8.6 °C, respectively (Chai et al., 2019).

Results

AS Event Identification and Analysis

According to description by Reddy (2007), alternative splicing events were sorted into 5 classes: Alternative 3′ splice site, Alternative 5′ splice site, Exon Skipping, Intron Retention and Mutually Exclusive Exons. The six samples grown in Xindu under normal conditions displayed an average of 205,496 AS events (Table 1; Table S5). The proportions of each AS type were analyzed in both zws-ms and zws-217. The two lines shared highly similar ratios of each AS type, with the alternative 5′ splice site and alternative 3′ splice site types being the most commonly observed, at 43.48% and 42.77% of AS events for both lines, respectively. The mutually exclusive exons type was the next most common (6.61%), followed by the Intron Retention type (5.92%), and the least common types was Exon Skipping, which represented just 1.22% of the AS events (Fig. 2A; Table 1).

Table 1:

Numbers of alternative splicing events in six samples.

Area	Line	Sample ID	Alternative 3′ splice site	Alternative 5′ splice site	Exon skipping	Intron retention	Mutually exclusive exons
Xindu	zws-ms	T01	87699	89187	2333	11981	12765
		T02	87603	88889	2404	11201	13100
		T03	88055	89623	2638	12184	14060
	zws-217	T04	87953	89443	2567	12501	13887
		T05	88012	89588	2581	12722	13901
		T06	88001	89408	2549	12373	13768
Ma’erkang(colder area)	zws-ms	T01	48432	49501	1970	14764	5299
		T02	46307	47031	1526	9560	3696
		T03	46594	47558	1712	12686	4464
	zws-217	T04	48204	49168	1596	12722	4733
		T05	46902	47912	1754	12194	4777
		T06	47759	48853	1822	13718	5105

DOI: 10.7717/peerj.10135/table-1

Notes:

T01, T02, and T03: Buds of three independent zws-ms plants at the budding stage; T04, T05, and T06: Buds of three independent zws-217 plants at the budding stage.

Alternative 3′ splice site: different-size mRNAs are produced depending on the usage of a proximal or distal 3′ splice site; Alternative 5′ splice site: different-size mRNAs are produced depending on the use of a proximal or distal 5′ splice site; Exon Skipping: an exon is either included or excluded from the mRNA; Intron Retention: an intron is either retained or excised in the mRNA, resulting in different-size transcripts; Mutually Exclusive Exons: adjacent exons are spliced in such a way that only one of them is included at a time in the mRNA.

As to the plants grown in colder area Ma’erkang, the two lines also shared highly similar ratios of each AS type: the alternative 5′splice site and alternative 3′splice site types represented the greatest proportion, at 42.13% and 41.29%, respectively; while the exon skipping accounted for the least proportion 1.51% (Fig. 2B; Table 1). The number of each type of AS events were significantly reduced in colder area Ma’erkang, except for Intron Retention.

Annotation of the alternatively spliced genes

To study the biological functions of the genes with AS events, GO and KEGG pathway enrichment analyses were performed. The GO annotations for AS genes from plants in Xindu included 17 terms involved in biological processes (BP; Fig. 3A), 17 terms associated with cellular components (CC), and 20 terms involved in molecular functions (MF). The most highly enriched BP terms observed in the alternatively spliced genes included “cellular process”, “single-organism process” and “metabolic process”. The most common CC categories were “cell”, “cell part” and “organelle”. In the MF category, the most enriched terms were “binding”, “catalytic activity” and “nucleic acid binding transcription factor activity”. Plants grown in Ma’erkang showed highly similar GO data to that in Xindu: 22 terms involved in BP (Fig. 3B), 15 terms associated with CC, and 15 terms involved in MF. Moreover, data from Ma’erkang and Xindu showed the same top-3 most enriched terms in each category.

These KEGG pathways were classified into five major groups: metabolism, genetic information processing, cellular processes, environmental information processing, and organismal systems. Of these, the subgroups “biosynthesis of amino acids”, “carbon metabolism”, “ribosome”, and “RNA transport” contained the highest number of annotated genes (Fig. 4A). Data from Ma’erkang showed similar subgroups containing the most of AS genes: “ribosome”, “carbon metabolism”, “biosynthesis of amino acids” and “plant hormones signal transduction” contained the highest number of annotated genes (Fig. 4B).

Discussion

As an important post-transcriptional metabolic event, AS is involved in many plant growth and developmental processes, such as flowering induction (Eckardt, 2002; Slotte et al., 2009) and the responses to environmental fluctuations and pathogen attacks (Barbazuk, Fu & McGinnis, 2008). To the best of our knowledge, AS events have seldom been reported to regulate the development of flower/fruit morphology in higher plants. This study is the first to analyze the role of AS events in rapeseed flower/fruit development as a whole, let alone those related to the multi-silique trait.

We previously described the morphology and inheritance of the multi-silique trait in B. napus (Jiang et al., 1998), investigating the associated regions of chromosomes at the genomic level and transcriptomically exploring the DEGs in multi-silique and single-silique plants (Chai et al., 2019). The multi-silique trait was found to be controlled by three recessive alleles and was significantly affected by environment; however, the mechanisms by which environmental factors affect this trait remained unknown, even if we knew that temperature could switch on/off the multi-silique trait (Chai et al., 2019). As mentioned above, AS is a pathway by which the environment could regulate plant physiology, therefore in this study, we analyzed AS events in order to investigate the mechanism by which plants perceive temperature fluctuations.

In this study, we sampled the buds of three individual plants from zws-ms and zws-217 lines in both Xindu and colder area Ma’erkang, and then subjected them to RNA-seq. All of the four groups generated sufficient data, which was validated by qPCR in earlier and present study successively. The samples in Ma’erkang group generated less data, mainly due to a reduced sequencing depth; nevertheless, this still provided enough data of high quality and assured the accuracy of the subsequent analysis.

We identified all of the genes with AS events in the zws-ms and zws-217 plants. In Xindu group, 11 AS genes were significantly differentially expressed between the multi-silique zws-ms line and its NIL, zws-217, which produces normal siliques. We analyzed their annotations and orthologs in Arabidopsis. One such ortholog, AT5G15470 (also known as Galacturonosyltransferase 14, GAUT14), is involved in cell wall pectin biosynthesis (Caffall et al., 2009), and the gaut13 gaut14 double mutant was previously shown to be defective in pollen tube growth (Wang et al., 2013a; Wang et al., 2013b). AT3G15420 (the ortholog of BnaA02g03080D) and AT3G10070 (the ortholog of BnaA09g45000D) encode subunits of the transcription factor complexes TFIIIC and TAF12, respectively. The former does not appear to be substantially involved in plant development; however, some members of the TAF family are involved in the regulation of morphology. The transgenic expression of TAF10 from clustered yellowtops (Flaveria trinervia) in Arabidopsis limited the development of the indeterminate inflorescence and resulted in the production of deformed leaves (Furumoto et al., 2005). By contrast, the taf mutant in Arabidopsis has abnormal phyllotaxis and lacks proper vegetative meristem activity (Tamada et al., 2007), indicating the important roles played by the TAFs in plant morphological development. Another DEG AS gene, BnaA04g16220D, is not annotated, and its Arabidopsis ortholog AT1G14800 is simply listed as an uncategorized nucleic acid-binding, OB-fold-like protein. The AS gene orthologs AT2G04900 and AT1G15060 encode an unknown protein and an uncategorized alpha/beta hydrolase family protein, respectively, so their roles in the regulation of the multi-silique trait are also currently unclear.

Another ortholog for differentially expressed AS gene, AT3G54620, is reported to encode a bZIP transcription factor-like protein. Members of this protein family are typically reported to regulate plant tolerance of environmental stresses. The transgenic expression of the maize (Zea mays) gene ZmbZIP72 in Arabidopsis enhanced its drought and salt tolerance (Ying et al., 2012), while BnbZIP3, a ramie (Boehmeria nivea) bZIP transcription factor, also increased the drought, salinity, and heavy metal tolerances of transgenic Arabidopsis (Huang et al., 2016). These genes are also involved in the regulation of other processes; for example, the repression of a bZIP transcription factor gene OsABI5 expression in rice resulted in low fertility (Zou et al., 2008), while the transgenic expression of tomato (Solanum lycopersicum) SlbZIP2 in tobacco (Nicotiana benthamiana) increased leaf thickness (Seong et al., 2016). To date, however, there are no reports of bZIP genes playing a significant role in flower/fruit morphology.

Other AS gene orthologs included AT5G16210, encoding a member of the HEAT repeat-containing protein family, which are considered to be involved in intracellular transport (Hernández-Torres, Jaimes-Becerra & Chomilier, 2014; Oeffinger, Dlakic & Tollervey, 2004). Although BnaC06g16950D is not annotated, its ortholog, AT3G59000, was identified as encoding an F-box/RNI-like superfamily protein in Arabidopsis, which typically function in the plant hormone signaling pathways (Gao et al., 2009). Similarly, the ortholog AT4G16900 encodes a TIR-NBS-LRR class protein, which are known to be involved in disease resistance (Xun et al., 2019) and hormonal responses (Sarazin et al., 2015). Moreover, Hewezi et al. (2006) unexpectedly found that these proteins are associated with developmental abnormalities; transgenic sunflowers (Helianthus annuus) expressing the antisense sequence complementing PLFOR48, which encodes a TIR-NBS-LRR-type protein, showed stunted growth and a reduction in apical dominance; whereas the pods of transgenic tobacco (N. tabacum) lacking PLFOR48 expression were smaller and showed severe deformations. This indicates that TIR-NBS-LRR-type proteins can regulate the morphology of plants, including fruit morphology, to some extent. Finally, AT1G10760, the ortholog of AS gene BnaC08g49610D, which encodes a GWD protein required for starch degradation, is involved in carbohydrate metabolism (Nadolska-Orczyk et al., 2017). This gene was also reported to regulate seed size; Pirone et al. (2017) found that the length and width of the mature seeds were reduced in the gwd1 Arabidopsis mutant, while their density was increased.

To summarize, AT5G15470, AT3G10070, AT3G54620, AT4G16900, and AT1G10760 are all known to be involved in plant development; therefore, their corresponding rapeseed orthologs, BnaA02g02630D, BnaA09g45000D, BnaAnng30260D, BnaC07g33980D, and BnaC08g49610D, the expression levels of which differed significantly between zws-ms and zws-217, are considered to be potential candidate genes regulating the multi-silique trait.

After that, we continued to investigate data from Ma’erkang, where it is colder and the multi-silique trait in zws-ms line disappeared. Due to the importance of the above-mention 11 differentially expressed AS genes, we first paid attention to them and found that three of them had the same expression models as in Xindu. In other words, they were independent of temperature; in addition, combined with their annotations, they were excluded from the potential candidate genes responding to environmental factors. On the other hand, other 8 AS genes stopped being differentially expressed between zms-ms and zws-217 in Ma’erkang. That is to say, the expression level of these 8 AS genes (BnaA02g02630D, BnaA02g03080D, BnaA04g16220D, BnaA09g45000D, BnaA09g45260D, BnaC02g06440D, BnaC07g33980D and BnaC08g49610D) were environment-specific. Besides, BnaC08g36010D, BnaC08g39130D and BnaC08g39360D, which were annotated to flower/ovule/fruit-related terms, these were differentially expressed between zws-ms and zws-217 specifically in Ma’erkang. Moreover, we also found nine temperature-responding AS genes differentially expressed in the colder area: BnaAnng17190D, BnaC01g27600D, BnaC08g39130D, BnaC09g53990D, BnaA07g19340D, BnaA05g28590D, BnaC06g15710D, BnaC08g36010D and BnaCnng24040D. The lower temperature motivated more responding genes and this also explained the reason for increased number of AS genes clustered in each KEGG pathway.

We also explored the line-specific AS genes, which were similarly expressed between zws-ms and zws-217, but contained stable and particular AS event(s) that differed between these two lines. These genes are likely to qualitatively regulate the multi-silique trait. In this case, we could obtain better results by fine-classification of the AS types into 12 subclasses, rather than the five classes mentioned above. Fine classifications could better identify differences between AS types more precisely and subtly. Thus, we found 205 genes of this type, of which 187 could be annotated. Due to the rarity of the multi-silique trait, we did not obtain much useful information from the KEGG pathway analysis. This meant that we were unable to relate this metabolic pathway information to the multi-silique trait directly; however, the GO analysis provided more potential clues. Among these, eight genes were considered to be associated with flower/carpel/ovule development. BnaC06g32640D is annotated as being involved in the regulation of the vegetative-to-reproductive phase transition in the meristem (GO:0010228) and in ovule development (GO:0048481). The IR of it in zws-217 seemed to block the multi-silique, while the TSS occurred in zws-ms as if it was positively related to multi-silique trait. Its Arabidopsis ortholog, AT1G71692, is annotated as AGAMOUS-LIKE12 (AGL12). Peng et al. (2015) isolated the BnFUL gene in rapeseed, which is homologous to AGL8 in Arabidopsis. Although BnFUL was hypothesized to be involved in enhancing pod-shattering resistance, when introduced into Arabidopsis, two of the five transgenic plants expressing BnFUL unexpectedly had a multi-silique phenotype. However, the mechanisms by which BnFUL generates this multi-silique phenotype remain elusive thus far, making the AGL12 gene identified in this study a potentially important candidate gene. The ortholog gene of BnaCnng68400D, AT5G15020, encodes an SIN3-LIKE 2 protein (SNL2) known to be important for seed germination or dormancy (Wang et al., 2016; Wang et al., 2013a; Wang et al., 2013b). The zws-ms line-specific AS events of it were detected in Xindu, where the two NILs were distinct from each other in the multi-silique trait, while in Ma’erkang, where the distinguishing trait was switched off, these AS event disappeared in zws-ms. This showed a positive correlation to the trait.

The other six selected line-specific AS genes showed particular AS events in zws-217 in Xindu, implying potential inhibition of the multi-silique morphology. However, these AS events did not occur in zws-ms or zws-217 in Ma’erkang, representing unknown complexity of the mechanisms, which were strongly influenced by environment. Fortunately, their orthologs in Arabidopsis provided many useful clues. AT2G20180 and AT4G00050, orthologs of BnaC07g00780D and BnaC03g32190D respectively, both encode phytochrome interacting factors (PIFs). Several transcription factors (AP1, SVP, LFY, AG, and SEP3) involved in the regulation of flowering are known to bind to the PIFs, suggesting a direct link with the reported flowering phenotype of the pif mutants (Leivar & Monte, 2014). AT5G17270 encodes a prenylyltransferase superfamily protein; however, to the best of our knowledge, there have been no reports about its development-related functions. The Arabidopsis ortholog of BnaC05g34570D is AT3G18600, which encodes a P-loop-containing nucleoside triphosphate hydrolase. While few studies have reported the functions of these proteins, Liu et al. (2016) reported that, in sesame (Sesamum indicum), one gene encoding a P-loop-containing nucleoside triphosphate hydrolase showed a reduced expression level in sterile buds, indicating that they may play a role in specifying/determining tapetal fate and development. Another line-specific AS ortholog, AT3G54660, encodes a glutathione reductase (GR), which was found to increase the fineness (mass per unit length) and bundle strength of cotton (Gossypium hirsutum) fiber when transgenically expressed (Tuttle et al., 2015). Since cotton fibers are single cells initiating from the epidermis of the outer integument of the ovules (Ruan et al., 2004), it can be inferred that GR regulates ovule development to some extent. AT3G28730, also known as structure-specific recognition protein SSRP1, was also found to regulate floral development, as the ssrp1-2 mutant Arabidopsis produced small and deformed petals with shorter stamens (Lolas et al., 2010).

To date, there is some evidence to show that the line-specific AS orthologs AT2G20180, AT4G00050 and AT3G54660 are related to the regulation of flower/fruit morphology, with clear roles reported for AT1G71692 and AT3G28730. Consequently, their orthologs in rapeseed, BnaC06g32640D and BnaC07g25280D, respectively, are considered to be important candidate genes regulating the multi-silique trait by conferring or removing some specific line-specific AS events in varied environments. In addition, BnaCnng68400D, of which AS events represented a positive correlation to morphology with or without multi-siliques, was also noteworthy.

Some of the genes/loci controlling silique development in Brassica plants have previously been reported. In addition to those regulating traits such as the seed weight and silique length (Liu et al., 2015) and the number of seeds per silique in B. napus (Li et al., 2015), some genes related to silique morphology have been cloned and functionally analyzed. Xiao et al. (2013) fine-mapped a multi-locular silique gene, Bjln1, to a 208-kb region on chromosome A7 in Brassica juncea and then revealed that it was the mutations in the CDS and promoter of BjuA07.CLV1 gene (equivalent to Bjln1) to cause the multi-locular trait (Xiao et al., 2018). Both Fan et al. (2014) and Yadava et al. (2014) reported that a mutation in BrCLV3, a homologue of CLAVATA3 in Arabidopsis, caused the production of multi-locular siliques in B. rapa. However, the multi-silique (or multi-pistil) phenotype of zws-ms is different from the above-motioned multi-locular trait; zws-ms produces three pods on each carpopodium, rather than multiple loculi per pod.

Few studies have investigated this multi-silique trait in rapeseed; however, there have been similar reports of multi-pistil traits in other crops, particularly in wheat (Triticum aestivum): Duan et al. (2015) discovered a male-sterile wheat mutant, dms, with a dwarf status and multi-pistils, a pleiotropic phenotype found to be controlled by a single recessive gene, which was not identified. Guo et al. (2019a) and Guo et al. (2019b) reported another multi-ovary trait in the wheat line DUOII, which was controlled by a dominant gene, and used a proteomics approach to propose some candidate proteins. Yang et al. (2017) mapped a gene promoting the formation of three pistils (Pis1) to chromosome 2D and identified some candidate genes according to their annotations, while Zhu et al. (2019) discovered a wheat multi-pistil mutant, 12TP, which was found to contain a semidominant mutation located on chromosome arm 2DL. Although several studies have explored the multi-pistil trait in wheat, no one has identified any of the specific genes responsible yet.

To sum up, the seven candidate genes mentioned above, including the five differentially expressed AS genes of interest (BnaA02g02630D, BnaA09g45000D, BnaAnng30260D, BnaC07g33980D, and BnaC08g49610D) and the two genes with line-specific AS events (BnaC06g32640D and BnaC07g25280D), are therefore hypothesized to regulate the multi-silique trait in rapeseed zws-ms, based on their AS expression levels or line-specific AS events altered by environment. These findings lay a foundation for further functional analyses for future experiments.

Conclusions

The utilization of heterosis is a way to increase the yield or improve the quality of crops. Exploring new germplasm resources and genes, as well as clarifying their inheritance, is the foundation of obtaining excellent hybrids. This study provides a novel inspection into the multi-silique trait in rapeseed from the transcriptional perspective by AS responding to environment, deepening the understanding of its molecular mechanism. Further functional verifications are now undergoing study.

Supplemental Information