Differential Alternative Splicing Genes in Response to Boron Deficiency in Brassica napus

Alternative splicing (AS) can increase transcriptome diversity, protein diversity and protein yield, and is an important mechanism to regulate plant responses to stress. Oilseed rape (Brassica napus L.), one of the main oil crops in China, shows higher sensitivity to boron (B) deficiency than other species. Here, we demonstrated AS changes that largely increased the diversity of the mRNA expressed in response to B deficiency in B. napus. Each gene had two or more transcripts on average. A total of 33.3% genes in both Qingyou10 (QY10, B-efficient cultivar) and Westar10 (W10, B-inefficient cultivar) showed AS in both B conditions. The types of AS events were mainly intron retention, 3′ alternative splice site, 5′ alternative splice site and exon skipping. The tolerance ability of QY10 was higher than that of W10, possibly because there were far more differential alternative splicing (DAS) genes identified in QY10 at low B conditions than in W10. The number of genes with both DAS and differentially expressed (DE) was far lower than that of the genes that were either with DAS or DE in QY10 and W10, suggesting that the DAS and DE genes were independent. Four Serine/Arginine-rich (SR) splicing factors, BnaC06g14780D, BnaA01g14750D, BnaA06g15930D and BnaC01g41640D, underwent differentially alternative splicing in both cultivars. There existed gene–gene interactions between BnaC06g14780D and the genes associated with the function of B in oilseed rape at low B supply. This suggests that oilseed rape could regulate the alterative pre-mRNA splicing of SR protein related genes to increase the plant tolerance to B deficiency.


Introduction
Alternative splicing (AS) is a widespread phenomenon in the majority of eukaryotic organisms that can generate one or multiple mRNAs isoforms from the same precursor mRNA (pre-mRNA) by using different splice sites [1,2]. The basic splicing process includes assembly of spliceosome and splicing [3]. Splice sites recognition of constitutive splicing by spliceosomes is regulated by 5 and 3 consensus cis-sequences and branchpoint [1]. Nevertheless, alternative splicing splice-site selection is also determined by cis-acting elements and trans-acting factors, which contain exonic/intronic

Plant Materials and Growth Conditions
B-efficient B. napus cultivar Qingyou 10 (QY10) and B-inefficient B. napus cultivar Westar 10 (W10) were used in this study. The seedlings were cultured using Hoagland solution [30] in a growth chamber. Then, 25 µL H 3 BO 3 (10 µmol/L) and 0.25 µL H 3 BO 3 (10 µmol/L) were applied in the solution to create a treatment with sufficient B and a treatment with deficient B with three replicates, respectively. The seedlings were grown in an illuminated culture room at 20 • C under a 16 h light/8 h dark cycle. The plants were arranged on a plant tray that had 9 × 6 holes with 5.6 cm between plants. The nutrient solution was oxygenated by an air pump and was replaced weekly. The pH of the nutrient solution was maintained at 5.7-6.0. When the seedlings showed B deficient symptoms (30 days), the root (R), old leaves (OL) and juvenile leaves (JL) of three plants of QY10 and W10 were sampled for RNA extraction, respectively. All samples were immediately frozen in liquid nitrogen and then stored at −80 • C.

RNA-Seq
Total RNA was extracted using an RNeasy Plant Kit (Bioteke RP1202) according to the manufacturer's instructions. First, 1.5 µg total RNA from each sample were used for preparing RNA-Seq libraries (cDNA libraries). The manufacturer's protocol (the Illumina Truseq RNA sample prep Kit, Illumina Inc., San Diego, CA, USA) was used to generate sequencing libraries. The PCR amplification products were checked, excised and purified by agarose gel electrophoresis and MinElute PCR Purification Kit (QIAGEN, Hilden, Germany). The final products were loaded onto flow cell channels with a concentration of 2 pM by using 2 × 100 bp pair-end sequencing strategy. The sequencing analysis were conducted in the National Key Laboratory of Crop Genetic Improvement, Huazhong Agricultural University using Illumina Hiseq TM 2000 platform (Illumina Inc., San Diego, CA, USA). All raw data were deposited in the NCBI Sequence Read Archive (Bioproject accession number, PRJNA393069) [31].

Reads Alignment, Transcript Assembly and Junction Prediction
The average reads generated from these RNA-Seq libraries was more than 170 million (Table S1). All the above sequences were trimmed based on quality (Q ≥ 30) and 68.2% of clean reads were mapped on the reference genome of Darmor-bzh (Table 1) using TopHat software (v2.0.12) [32]. The aligned reads were assembled into transcripts using cufflinks software (v2.2.1) [33]. The cuffcompare package with cufflinks suite program was used to identify novel genes and transcripts that are identical to the B. napus reference genome. In addition, TopHat software (v2.0.12) was also used to predict the splice junctions, and to identify the known and novel splice junctions according to the gene information with default parameters. In principal component analysis (PCA) of all genes, the replicate samples showed a high similarity with respect to the first two principal components, little variation within group and a good separation of groups ( Figure S1). Table 1. The starting reads, clean reads ratio, and mapped reads ratio against in the B. napus reference genome of c.v. Darmor-bzh.

Detection of Alternative Splicing Events
All RNA CEL files were normalized using the RPKM (Reads Per kb per Million reads) method [34]. Alternative splicing events were identified by the AS transcriptional landscape visualization tool (AStalavista) [35]. Source codes for the AStalavista software (v3.0) are available online at http:// genome.crg.es/astalavista/. TopHat (v2.0.12) was used to calculate the AS events numbers. This study focused on four main types of AS: intron retention, exon skipping (cassette exons), alternative 3 splice site (alternative acceptor) and alternative 5 splicing site (alternative donor).

Identification of Differential Alternative Splicing Events and Differentially Expressed Genes
Differential alternative splicing (DAS) events were analyzed using the program Replicate Multivariate Analysis of Transcript Splicing (rMATS, v3.0.9) [36]. A T-test was used to calculate p value of splicing events with ∆ ψ. A stringent threshold, p ≤ 0.05 and |∆ ψ| ≥ 0.05, was used to define the DAS events of AS events between different cultivars or different tissues of oilseed rape. The events with p values less than 0.05 were identified as significantly different events. Differentially expressed genes (DEGs) were calculated using a t-test (corrected by Benjamin-Hochberg false discovery rate (FDR) multiple testing). Probe-sets with a FDR corrected p-value ≤ 0.05 and fold change of >2 were considered to be differentially expressed.

Semi-Quantitative RT-PCR Analysis
Total RNA was extracted using an RNeasy Plant Kit (Bioteke RP1202) according to the manufacturer's instructions. The DNAase treatment is mentioned above. Synthesizing the cDNA used a reaction solution containing 1 µg RNA with HiFiScript, dNTP Mix, Primer Mix, 5 × RT Buffer, DTT and HiFiScript. The solution was incubated at 42 • C for 50 min, followed by 85 • C for 5 min. The cDNA was used as a PCR template in a 20 µL reaction system of semi-quantitative RT-PCR. To validate the DAS events, 12 pairs of primers were designed based on different DAS events including intron retention, alternative 3 splice site and alternative 5 splice site. The primers used are listed in Table S2.

Gene Ontology (GO) Analysis and KEGG Analysis
Blast2GO was used to determine the gene functional category with a cut-off of 1E-5 [37]. An internal Perl script was used to perform GO annotation based on Open Biological and Biomedical Ontologies (OBO) (http://purl.obolibrary.org/obo). Biochemical pathway was analyzed by Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology based on annotation system version 2.0 (KOBAS) [38]. The different pathways shown in the figures were chosen based on statistical significance (p < 0.05).

Novel Transcripts and Novel Genes in B. napus
Assembly of QY10 and W10 under B sufficient and deficient conditions identified 135,036 expressed transcripts and 56,203 expressed genes, and each gene had two or more transcripts on average (Table 2). Many new transcripts and genes were discovered by comparing the transcripts of assembling and the reference genome. The numbers of novel isoforms were significantly higher than that of the reference isoforms in both QY10 and W10 whether under B sufficient or deficient conditions ( Figure 1). GO analysis of the novel genes revealed that the categories of "ATP hydrolysis coupled proton transport" in biological process, "proton-transporting ATPase activity" in molecular function and "proton-transporting V-type ATPase, V0 domain" in cellular component were significantly enriched for common novel genes in QY10 and W10 under both B conditions ( Figure 2). Change of ATP activity could affect the stability of cell membranes by changing the proton gradient on the plasma membrane at low B.

Identification of AS Events in QY10 and W10 under B Sufficient and Deficient Conditions
In total, 11,040,252 splice junctions were identified in all 36 RNA-Seq libraries (Table S3). Four major types of AS events, namely exon skipping, alternative 5′ splice site, alternative 3′ splice site and

Identification of AS Events in QY10 and W10 under B Sufficient and Deficient Conditions
In total, 11,040,252 splice junctions were identified in all 36 RNA-Seq libraries (Table S3). Four major types of AS events, namely exon skipping, alternative 5′ splice site, alternative 3′ splice site and

Identification of AS Events in QY10 and W10 under B Sufficient and Deficient Conditions
In total, 11,040,252 splice junctions were identified in all 36 RNA-Seq libraries (Table S3). Four major types of AS events, namely exon skipping, alternative 5 splice site, alternative 3 splice site and intron retention, were investigated in this study ( Figure 3). Under B sufficient conditions, the numbers of AS events and genes in the root, juvenile leaves and old leaves of W10 were higher than those of QY10 ( Figure 4). The numbers of AS events and genes in root and juvenile leaves of W10 were also higher than those of QY10 under low B conditions ( Figure 4). Moreover, the numbers of AS events and genes in both QY10 and W10 under B deficient conditions were higher than those under B sufficient conditions. Genes 2019, 10, x FOR PEER REVIEW 7 of 22 intron retention, were investigated in this study ( Figure 3). Under B sufficient conditions, the numbers of AS events and genes in the root, juvenile leaves and old leaves of W10 were higher than those of QY10 ( Figure 4). The numbers of AS events and genes in root and juvenile leaves of W10 were also higher than those of QY10 under low B conditions ( Figure 4). Moreover, the numbers of AS events and genes in both QY10 and W10 under B deficient conditions were higher than those under B sufficient conditions.  The AS genes with the types of splicing of "intron retentions (RI)" and "alternative 3′ splice site (A3SS)" were the largest fraction, whereas the AS genes with "exon skipping (ES)" were the least fraction among all the AS genes in both cultivars whether under B sufficient or B deficient conditions ( Figure 5). The number of AS genes with "alternative 5′ splice site (A5SS)" was between the types of "RI", "A3SS" and "ES". Most AS genes had one type of AS in the root, juvenile leaves and old leaves intron retention, were investigated in this study ( Figure 3). Under B sufficient conditions, the numbers of AS events and genes in the root, juvenile leaves and old leaves of W10 were higher than those of QY10 ( Figure 4). The numbers of AS events and genes in root and juvenile leaves of W10 were also higher than those of QY10 under low B conditions ( Figure 4). Moreover, the numbers of AS events and genes in both QY10 and W10 under B deficient conditions were higher than those under B sufficient conditions.  The AS genes with the types of splicing of "intron retentions (RI)" and "alternative 3′ splice site (A3SS)" were the largest fraction, whereas the AS genes with "exon skipping (ES)" were the least fraction among all the AS genes in both cultivars whether under B sufficient or B deficient conditions ( Figure 5). The number of AS genes with "alternative 5′ splice site (A5SS)" was between the types of "RI", "A3SS" and "ES". Most AS genes had one type of AS in the root, juvenile leaves and old leaves The AS genes with the types of splicing of "intron retentions (RI)" and "alternative 3 splice site (A3SS)" were the largest fraction, whereas the AS genes with "exon skipping (ES)" were the least fraction among all the AS genes in both cultivars whether under B sufficient or B deficient conditions ( Figure 5). The number of AS genes with "alternative 5 splice site (A5SS)" was between the types of "RI", "A3SS" and "ES". Most AS genes had one type of AS in the root, juvenile leaves and old leaves in W10 and QY10 under B sufficient and deficient conditions ( Figure 6 and Table S4). Only a small portion of AS genes had all four types of AS in the same organ across two cultivars and two B treatments ( Figure 6 and Table S4). The number of AS genes with the four types of AS in the juvenile leaves of W10 under B deficient condition was the most and that in the juvenile leaves of QY10 under B sufficient condition was the least. The number of AS genes with both "RI" and "A3SS" was the greatest and that of AS genes with both "ES" and "A5SS" was the least in each part investigated ( Figure 6 and Table S4). treatments ( Figure 6 and Table S4). The number of AS genes with the four types of AS in the juvenile leaves of W10 under B deficient condition was the most and that in the juvenile leaves of QY10 under B sufficient condition was the least. The number of AS genes with both "RI" and "A3SS" was the greatest and that of AS genes with both "ES" and "A5SS" was the least in each part investigated ( Figure 6 and Table S4).
The categories of "cell wall modification", "glucose-6-phosphate transport", "ion transport", and "response to stress" were mostly enriched for AS genes, while the categories of "DNA packaging", "ribosomal export" and "protein location" were enriched for no AS genes ( Figure 7). These demonstrated that the genes with AS and non-AS had different functions. Moreover, the AS genes of QY10 and W10 identified in B deficient conditions showed more than that in B sufficient conditions. The AS of most of the genes tended to occur under abiotic stress conditions (Figure 7).   The categories of "cell wall modification", "glucose-6-phosphate transport", "ion transport", and "response to stress" were mostly enriched for AS genes, while the categories of "DNA packaging", "ribosomal export" and "protein location" were enriched for no AS genes (Figure 7). These demonstrated that the genes with AS and non-AS had different functions. Moreover, the AS genes of QY10 and W10 identified in B deficient conditions showed more than that in B sufficient conditions. The AS of most of the genes tended to occur under abiotic stress conditions (Figure 7).

Differential Alternative Splicing in QY10 and W10 in B Deficient Conditions
In total, 159, 190 and 163 DAS genes were identified in root, juvenile leaves and old leaves of QY10, respectively, under B deficient conditions ( Figure 8). Among them, 33 DAS genes were detected simultaneously in root, juvenile leaves and old leaves, and were involved in several important biological processes, such as "Phosphoglucomutase", "Translation elongation factor EF1B", "glutathione S-transferase phi 8" and "calcineurin B-like protein 9" (Table 3). In total, 24, 40 and 74 DAS genes were distinguished in root, juvenile leaves and old leaves in W10, respectively, in B deficient condition. Only two DAS genes coexisted in root, juvenile leaves and old leaves ( Figure  8). They functioned as cystatin and phosphoribulokinase, respectively (Table 3). A total of 63 DAS genes were detected simultaneously in all three organs of QY10 and W10 under B deficient conditions. Moreover, 313 and 58 cultivar-specific DAS genes were identified in QY10 and W10, respectively ( Figure 8).
The percentages of the four major types of DAS genes in root, juvenile leaves and old leaves of QY10 under B deficiency, and juvenile leaves and old leaves of W10 under B deficiency, all showed

Differential Alternative Splicing in QY10 and W10 in B Deficient Conditions
In total, 159, 190 and 163 DAS genes were identified in root, juvenile leaves and old leaves of QY10, respectively, under B deficient conditions (Figure 8). Among them, 33 DAS genes were detected simultaneously in root, juvenile leaves and old leaves, and were involved in several important biological processes, such as "Phosphoglucomutase", "Translation elongation factor EF1B", "glutathione S-transferase phi 8" and "calcineurin B-like protein 9" (Table 3). In total, 24, 40 and 74 DAS genes were distinguished in root, juvenile leaves and old leaves in W10, respectively, in B deficient condition. Only two DAS genes coexisted in root, juvenile leaves and old leaves (Figure 8). They functioned as cystatin and phosphoribulokinase, respectively (Table 3). A total of 63 DAS genes were detected simultaneously in all three organs of QY10 and W10 under B deficient conditions. Moreover, 313 and 58 cultivar-specific DAS genes were identified in QY10 and W10, respectively ( Figure 8).
that RI > A3SS, A5SS > ES. However, in the root of W10, the percentage of "ES" DAS genes was the highest, and that of "A5SS" DAS genes was the least among the four major types of DAS genes ( Figure 9). If we used W10 as control, the total percentages of "RI" and "A3SS" DAS genes were the highest, whereas that of "ES" DAS genes was the least in QY10 under B sufficient or deficient conditions.  Table 3. Thirty-three and two differential alternative splicing genes detected in root, juvenile leaves and old leaves simultaneously in B. napus cultivars Qingyou10 and Westar10, respectively.

Cultivars
GeneID Gene Description DAS-Type  Phosphoglucomutase/phosphomannomutase family protein RI Ribosomal L29 family protein ES BnaC09g54460D Ribosomal protein S13/S18 family A3SS BnaA07g25750D RNA-binding (RRM/RBD/RNP motifs) family protein RI TLD-domain containing nucleolar protein RI BnaA03g07610D Translation elongation factor EF1B A3SS BnaC04g56630D unknown protein RI BnaCnng63660D unknown protein RI BnaC03g03780D VND-interacting 1 (VNI1) A3SS The percentages of the four major types of DAS genes in root, juvenile leaves and old leaves of QY10 under B deficiency, and juvenile leaves and old leaves of W10 under B deficiency, all showed that RI > A3SS, A5SS > ES. However, in the root of W10, the percentage of "ES" DAS genes was the highest, and that of "A5SS" DAS genes was the least among the four major types of DAS genes ( Figure 9). If we used W10 as control, the total percentages of "RI" and "A3SS" DAS genes were the highest, whereas that of "ES" DAS genes was the least in QY10 under B sufficient or deficient conditions. BnaCnng63660D unknown protein RI BnaC03g03780D VND-interacting 1 (VNI1) A3SS

W10
BnaCnng40950D Cystatin/monellin superfamily protein A3SS BnaC07g51220D nicotinate phosphoribosyltransferase 1 (NAPRT1) A5SS Note: DAS, Differential alternative splicing; R, root; JL, juvenile leaves; OL, old leaves; A3SS, Alternative 3′ splice site; A5SS, Alternative 5′ splice site; RI, retain intron; ES, exon skip.  KEGG pathways showed that these DAS genes were involved in all kinds of stress response pathways in B. napus, such as "Glycolysis/Gluconeogenesis", "Calcium signaling pathway", "MAPK signaling pathway" and "Peroxisome" (Table 4). Some important DAS genes were involved in the B deficient responses, especially in the root in both QY10 and W10. For example, BnaC05g18490D was associated with the Glycolysis process, which encoded the phosphoglucomutase and regulated cell wall synthesis. B deficiency led to the increase in retention of intron 13 in the juvenile leaves and old leaves in QY10, and the decrease in retention of intron 13 in the old leaves in W10. In addition, BnaC07g18010D and BnaA02g29830D were involved in the cell wall synthesis and integrity of cell membrane, which is an important biological process in response to B deficiency. The two genes showed decreased intron retention in old leaves in W10 in B deficient conditions. The results of qRT-PCR of these genes are consistent with the RNA-Seq data ( Figure 10). The change of pre-mRNA splicing could influence the metabolic processes in QY10 and W10 in response to B deficiency.

DAS Genes and DE Genes in QY10 and W10 in B Deficient Conditions
Combined analysis of DE genes and DAS genes in root, juvenile leaves and old leaves of QY10 and W10 indicated that about 0.06% of DE genes showed AS (Table 5). Under B deficient conditions, Figure 10. Validation of DAS events in the root, juvenile leaves and old leaves of B. napus cultivars Qingyou10 (QY10) and Westar10 (W10) under boron deficient conditions. The band with the red asterisk showed the fragment generated by the alternative splicing of the gene, which was determined according to the size of the retained intron fragment. The forward and reverse primers, P1 and P2, were designed based on the exon between the retained intron. Bs, boron sufficient condition; Bd, boron deficient condition; QR, root of QY10; QOL, old leaves of QY10; QJL, juvenile leaves of QY10; WR, root of W10; WOL, old leaves of W10; WJL, juvenile leaves of W10. RI, retain intron.
Most of the intron retention transcripts that had a premature termination codon (PTC) are degraded by the "nonsense-mediated mRNA decay (NMD)" surveillance pathway or are targeted by the microRNAs [12,14]. Sequence analysis of the intron retention transcripts of QY10 and W10 at low B showed that PTCs were found in all of these transcripts. The decrease or increase in retention of intron could increase or decrease the abundance of the functional transcripts, respectively. For example, the functional transcript levels of BnaC05g18490D, BnaC07g18010D and BnaA02g29830D were up-regulated in W10 under B deficient condition. In contrast, the functional transcript levels of these genes were down-regulated or unchanged in QY10 under B deplete condition.

DAS Genes and DE Genes in QY10 and W10 in B Deficient Conditions
Combined analysis of DE genes and DAS genes in root, juvenile leaves and old leaves of QY10 and W10 indicated that about 0.06% of DE genes showed AS (Table 5). Under B deficient conditions, the number of genes with both DAS and DE in the root of QY10 was the highest (32), and in the juvenile leaves of W10 was the least (13) ( Table 5). If W10 were used as a control, the number of genes with both DAS and DE in the old leaves of QY10 under B sufficient conditions was the highest (31). No gene showed both DAS and DE in the juvenile leaves and old leaves of QY10 under B deplete conditions ( Table 5). The number of genes with both DAS and DE was lower than that of the genes that were either with DAS or DE in QY10 and W10, which demonstrated that the DAS and DE genes were independent. Table 5. The number of differential alternative splicing genes and differential expressed genes in the root, juvenile leaf and old leaf of B. napus cultivars Qingyou10 and Westar10 under boron deficient conditions. Functional categorization of the DAS and DE genes in the root, juvenile leaves and old leaves of QY10 and W10 under B deficient conditions revealed that most of DAS and DE genes enriched in different functional pathways in biological process. For example, the categories of "mRNA processing", "Succinyl-CoA metabolic process", "signal transduction" and "cellular response to stress" were only enriched for DAS genes, whereas the categories of "cell wall formation", "cell wall pectin metabolic process" and "cell wall organization" were only enriched for DE genes (Figure 11). Only a small group of DAS and DE genes were annotated to the same function. For example, the categories of "starch metabolism process" and "glucose metabolism process" were enriched for both DAS genes and DE genes ( Figure 11).

Sample Names DAS Genes DE Genes Overlap
The majority of the DAS genes had differential intron retention (DIR), such as BnaA07g33860D (Sulphate anion transporter); and alternative 3 splice site (DA3SS), such as BnaA01g30340D (SANT/Myb domain) and BnaA06g17710D (alpha-glucan phosphorylase 2, PHS2) (Table S5). These genes also showed significant up-regulation or down-regulation under B deplete conditions (Table S5).

SR Splicing Factors in QY10 and W10 under B Deficient Condition
SR splicing factors play a crucial role to regulate pre-mRNA splicing in plants. A total of 34 SR splicing factors were identified in QY10 and W10 under B deficient conditions (Table S6). Four of these SR splicing factors, BnaC06g14780D, BnaA01g14750D, BnaA06g15930D and BnaC01g41640D, underwent DAS. Under B deplete conditions, A3SS occurred in the genes of BnaC06g14780D in the root and juvenile leaves; A5SS occurred in the genes of BnaA01g14750D in old leaves of QY10; and ES occurred in the gene of BnaC01g41640D in the root and old leaves of QY10 and BnaA06g15930D in the juvenile leaves of W10 (Table S6).
BnaC06g14780D showed differential A3SS in the second intron of the root and juvenile leaves in QY10 under different B conditions (Figure 12). In addition, its expression level under B deficient conditions was lower than that in B sufficient conditions. Compared with the B. napus reference

SR Splicing Factors in QY10 and W10 under B Deficient Condition
SR splicing factors play a crucial role to regulate pre-mRNA splicing in plants. A total of 34 SR splicing factors were identified in QY10 and W10 under B deficient conditions (Table S6). Four of these SR splicing factors, BnaC06g14780D, BnaA01g14750D, BnaA06g15930D and BnaC01g41640D, underwent DAS. Under B deplete conditions, A3SS occurred in the genes of BnaC06g14780D in the root and juvenile leaves; A5SS occurred in the genes of BnaA01g14750D in old leaves of QY10; and ES occurred in the gene of BnaC01g41640D in the root and old leaves of QY10 and BnaA06g15930D in the juvenile leaves of W10 (Table S6).
BnaC06g14780D showed differential A3SS in the second intron of the root and juvenile leaves in QY10 under different B conditions ( Figure 12). In addition, its expression level under B deficient conditions was lower than that in B sufficient conditions. Compared with the B. napus reference genome of Darmor-bzh, the A3SS imported a novel sequence. The exon insertion of BnaC06g14780D in the root and juvenile leaves in QY10 could increase the abundance of functional transcript in response to B deficiency ( Figure 12). BnaA01g14750D had lower A5SS in intron 2 in the root of QY10 under B deficient conditions than B sufficient conditions. A5SS also produced a novel sequence with PTC ( Figure 12). The exon insertion of BnaA01g14750D in the root of QY10 could also increase the abundance of the functional transcript under B deficient conditions. In addition, BnaA06g15930D and BnaC01g41640D were exon skipping, which constituted only a small portion of differentially expressed alternatively spliced genes in plants (data not shown). genome of Darmor-bzh, the A3SS imported a novel sequence. The exon insertion of BnaC06g14780D in the root and juvenile leaves in QY10 could increase the abundance of functional transcript in response to B deficiency ( Figure 12). BnaA01g14750D had lower A5SS in intron 2 in the root of QY10 under B deficient conditions than B sufficient conditions. A5SS also produced a novel sequence with PTC ( Figure 12). The exon insertion of BnaA01g14750D in the root of QY10 could also increase the abundance of the functional transcript under B deficient conditions. In addition, BnaA06g15930D and BnaC01g41640D were exon skipping, which constituted only a small portion of differentially expressed alternatively spliced genes in plants (data not shown).
BnaC06g14780D was a seed gene, and its downstream target genes were associated with the function of B [18], such as expansion protein (BnaA09g52970D) and bZIP transcription factor (BnaC04g52770D) ( Table 6 and Figure S2).  BnaC06g14780D was a seed gene, and its downstream target genes were associated with the function of B [18], such as expansion protein (BnaA09g52970D) and bZIP transcription factor (BnaC04g52770D) ( Table 6 and Figure S2).

Discussion
Alternative splicing, which generates multiple transcripts from the same gene, is an important modulator of gene expression that can increase proteome diversity and regulate mRNA levels [15]. RNA-Seq data revealed AS in 48% of genes in B. napus cultivar "Darmor-bzh" [39]. In this study, 30% and 35% of intron-containing genes underwent AS in QY10 and W10 under B sufficient and deficient conditions, respectively (Table S4). The number of AS gene identified in this study was lower than that in the cultivar "Darmor-bzh", possibly because the sequencing depth of the former is lower than the latter. The increase in AS genes of B. napus under B deficient conditions ( Figure 4) indicated that AS might be an important strategy of posttranscriptional regulation, and the increase in AS could improve the molecular plasticity of plants to adapt to abiotic stress. Recently, a higher proportion of genes are detected showing AS under salt stress in Arabidopsis [29], drought stress in maize [40], heat shock in the moss Physcomitrella patens [28] and high temperature in grape [11]. AS regulates the plant response to abiotic stress are largely by targeting the abscisisc acid (ABA) pathway [15]. However, the AS genes of QY10 and W10 identified in B deficient conditions were not associated with the ABA signaling ( Figure 7).
In this study, the number of AS genes with the types of "RI", "A3SS" and "A5SS" were much more than that with "ES" and others in B. napus, and RI was the most prominent types of AS ( Figure 5). In B. napus cultivar "Darmor-bzh", intron retention is also frequent (62%), whereas exon skipping is rare (3%) [39]. High proportions of AS genes with the types of "RI" are also found in maize [41] and Physcomitrella patens [42]. Although the intron retention is highly repressed by elevated temperature in Physcomitrella patens, the AS genes with the types of "RI" constituted the largest fraction of alternatively spliced genes [28]. Moreover, the genes responded to B deficient condition in this study, such as BnaC05g18490D, BnaC07g18010D and BnaA02g29830D, which were associated with cell wall synthesis and integrity of cell membrane, showed differential RI during B deficiency (Table S5). AS events based on different splicing types may lead to functionally relevant changes in the protein products [43]. For example, 4% of R2R3-MYB genes had undergone AS events in soybean, which generate a variety of transcripts to increase the complexity of transcriptome [44]. In this study, the categories of "cell wall modification", "glucose-6-phosphate transport", "ion transport", and "response to stress" were mostly enriched for AS genes, while the categories of "DNA packaging", "ribosomal export" and "protein location" were enriched for non-AS genes (Figure 7), which indicated that the genes with AS and non-AS play distinct physiological roles.
B. napus is highly susceptible to B deficiency [20]. There are significant genotypic differences in the response to low-B stress among different B. napus cultivars [21]. Under B deplete conditions, 159, 190 and 163 DAS genes were identified in the root, juvenile leaves and old leaves of QY10 (B-efficient cultivar), respectively; however, only 24, 40 and 74 DAS genes were identified in the root, juvenile leaves and old leaves of W10 (B-inefficient cultivar) (Figure 8). The increase of AS events occurred under abiotic stress conditions could enhance the tolerance ability of plants [15]. Further experiments should be conducted to confirm whether the tolerance ability of QY10 was higher than that of W10 belong to the DAS genes generated in B-efficient cultivar (QY10) were far more than that in B-inefficient cultivar (W10) under B deplete conditions.
The number of genes with both DAS and DE were far lower than that of the genes that were either DAS or DE in QY10 and W10 (Table 5), which demonstrated that the DAS and DE genes were independent. Functional categorization of DAS and DE genes of root, juvenile leaves and old leaves in QY10 and W10 under B deficient conditions revealed that most of DAS and DE genes enriched in differential functional pathways in biological process. However, only a small group of DAS and DE genes were annotated to the same function (Table 5, Figure 11). The majority of the DAS genes had differential intron retention (DIR) and alternative 3 splice site (DA3SS), which also showed significant up-regulation or down-regulation in B deplete condition (Table S5). These results demonstrate that both transcriptional regulation and posttranscriptional regulation contribute to B. napus adaption to B deficiency.
SR splicing factors play a crucial role in regulation of pre-mRNA splicing in plants [45]. In Arabidopsis, the AS pattern of members of SR splicing factors have been shown to change under various stress conditions, such as high light intensity, salinity and temperature stress [46][47][48]. This suggests that the stress-induced changes in SR splicing factors could in turn alter the factors of downstream targets in response to stress environments [48]. In this study, 34 SR splicing factors in QY10 and W10 were found under B deficient conditions (Table S6). Four of them, namely BnaC06g14780D, BnaA01g14750D, BnaA06g15930D and BnaC01g41640D, showed differential alternative spliced genes (Table S6). The decrease of novel sequence of BnaC06g14780D and BnaA01g14750D in the root of QY10 under B deplete conditions suggested that the increase of functional isoforms in the two SR genes might improve plant tolerance in response to B deficiency.

Conclusions
A total of~33.3% genes showed AS in B. napus cultivars QY10 and W10. The AS genes with the types of splicing of "intron retentions (RI)" and "alternative 3 splice site (A3SS)" in both cultivars showed the largest fraction whether under B deficient or sufficient conditions. Further experiments should be conducted to confirm whether the tolerance ability of QY10 was higher than that of W10 was attributed to far more DAS genes were identified in QY10 under low B conditions than in W10. To find the functional transcript of the genes, such as SR splicing factor BnaC06g14780D, which responds to low boron stress, would provide a new way to increase the ability of B. napus to cope with B deplete stress.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/3/224/s1, Figure S1: Principal component analysis on all the RNA-Seq. data of the root, old leaves and juvenile leaves of Qingyou10 (a) and Westar10 (b) under B sufficient and deficient conditions. Figure S2: The network of SR splicing factors and their target genes in B. napus under B deficient conditions. Table S1: Differential alternative splicing genes in response to boron deficiency in B. napus. Table S2: Primer sequences of actin and candidate genes in B. napus for Semi-Quantitative RT-PCR. Table S3: Splice junctions identified in the 36 RNA-Seq libraries of B. napus cultivars Qingyou10 (QY10) and Westar10 (W10). Table S4: The numbers of AS genes with different types of alternative splicing in B. napus cultivars Qingyou10 (QY10) and Westar10 (W10). Table S5: Annotation of the genes in transcriptional and posttranscriptional regulation in the root, juvenile leaves and old leaves of B. napus cultivars Qingyou10 and Westar10 under boron deficient conditions. Table S6: The types of differential alternative splicing of SR splicing factor in B. napus cultivars Qingyou10 and Westar10 under boron deficient conditions. Author Contributions: Conceptual and experiment designs were by L.S., W.L., C.L. and J.G.; Experiments were conducted by J.G. and W.L.; RNA-seq data analysis was performed by J.G.; Reagents/materials/analysis tools were contributed by L.S., J.R. and C.L., and the report was written by G.Y., G.D., S.W., F.X., A.C., C.L., X.Z. and L.S. All authors have commented on, read and approved the final manuscript.