Identification and characterization of BnDELLAs
We have identified 10 BnDELLA genes in B. napus using the known 5 A. thaliana DELLAs (GAI, RGA, RGL1, RGL2, RGL3) peptide sequences as queries and performed BLASTP search in the B. napus genome database (GENOSCOPE http://www.genos cope.cns.fr/brass icanapus/) . To further confirm the BnDELLA proteins integrity in the B. napus, we analyze the retrieved sequences in different B. napus cultivar genome browser (BnPIR, http://cbi.hzau.edu.cn/bnapus) and manually corrected the redundant sequence information of the BnDELLAs and named them according to their loci. Based on these methods, we found that each member in 5 AtDELLAs corresponds to multiple homologs in B. napus (Table 1). At the same time, 5 DELLA genes in B. rapa, 4 in B. oleracea, 9 in B. juncea, and 5 in B. nigra were identified by using the same methods. Furthermore, we found that 10 BnDELLAs members are derived from its progenitor B. rapa and B. oleracea. The chemical properties of the DELLA proteins in B. napus were identified by using the ExPASy (Table 1). The genomic sequence length of the BnDELLAs ranged from the 1524-1740bp, with a molecular weight varying from 55827.07 D to 63323.35 D. Moreover, the pI values of the BnDELLA proteins varied from 4.69 to 5.94, which shows that these proteins are highly acidic. Besides, all BnDELLA proteins showed a negative value of the grand average of hydrophobicity (GRAVY), indicating that BnDELLAs are hydrophilic. Subcellular localization signals of the 10 BnDELLA proteins were detected in the nucleus. The names of the BnDELLA genes and their locus id are also shown in (Table1).
Evolutionary relationship and gene structure analysis of BnDELLAs
The DELLAs evolutionary history among six Brassicaceae species A. thaliana, B. napus, B. rapa, B. oleracea, B. juncea and B. nigra were deduce using the neighbor-joining method. Based on the phylogenetic analysis, 38 DELLA genes in which 5 AtDELLAs, 10 BnDELLAs, 5 BrDELLAs, 4 BoDELLAs, 9 BjDELLAs, and 5 BniDELLAs were cluster into three groups according to the topologies and bootstrap support (Figure 1). Group I contain GAI and RGA clade, Group II holds RGL1 clade, Group III includes RGL2 and RGL3 clade. B. napusDELLA genes were relatively closer to the A. thaliana. However, B. napus and B. rapaDELLAs shows 100% similarity between each other. In addition, a homolog of AtRGL1 was not identified in the B. oleracea compared to those of B.napus, B.rapa, B. juncea, and B. nigra, which might be due to gene loss during the evolution process or maybe the emerging of genome gaps in B. oleracea. However, in B. napus, Group I, II, III had 4, 2, 4 DELLA members, respectively. DELLA genes grouped into the same subfamily are previously known to have distinct or overlapping functions [12, 22, 26, 34]. To recognize the DELLA genes family diversification in B. napus, we have implemented the GSDS web analysis by comparing the CDS and corresponding genomic sequences of AtDELLAs, BnDELLAs, BrDELLAs, BoDELLAs and BjDELLAs. As shown in (Figure 2), members of the DELLA genes among denoted species are highly conserved and intron-less with only one exon, which indicates BnDELLAs are the retrogenes. Moreover, the exon location of DELLAs among different phylogenetic related species is conserved, suggesting a similar evolutionary relationship. However, the length of exon among DELLA subfamily was different. For example, BnRGL1 exon length was smaller than other BnDELLAs members in Group I and Group III, indicating gene structure diversification. In summary, the gene structure of the DELLA genes from different Brassicaceae species are highly conserved with some difference in the exon length (Figure 2).
Multiple sequence alignment and analysis of BnDELLAs motifs
The putative sequences of the DELLA proteins from B. napus,A. thaliana, B. oleracea, B. rapa , B. juncea, and B. nigra were aligned to explore amino acid conservation in B. napus. Based on the alignment, we found five homologues DELLA proteins from A. thaliana show higher percent amino acid similarity with B. napus (Table S1). Similar to A. thaliana, the B. napus and other denoted species contain highly conserved DELLA and GRAS domain at the N-terminal and C-terminal region, respectively. The N-terminal DELLA domain is shown to be involved in stabilizing the DELLA genes activity [35, 36], while GRAS domain acts as coregulators to interact with several transcriptional factors and regulators [28, 37, 38] (Figure S1). Moreover, the BnDELLAs amino acid sequence’s alignment showed a highly conserved consensus sequence at the C-terminal region, which possibly suggests an evolutionary relationship between BnDELLA genes (Figure S1). In addition, the presence of the (VHIID-PFYRE-SAW) and two leucine heptad repeats LHRI and LHRII on the C-terminal of the GRAS domain are responsible for the protein interaction [39, 40]. However, some studies have also proposed DELLA domain lower-affinity with intrinsically unstructured proteins [37, 41]. Overall, the domain arrangements of the DELLA members in B. napus are equivalent to other species. In addition, the secondary structure feature (alpha-helix and Beta sheets) from the AtRGL1 accessions number (At1G66350.1) was displayed in (Figure S1). However, the predicted secondary structures of all DELLA genes from the denoted plant species were relatively different.
Furthermore, to gain more insights into the putative roles of BnDELLAs diversity in B. napus, we generated a graph showing domains and their position on AtDELLAs and BnDELLA protein members. We found that the DELLA and GRAS domains are conserved in all DELLA proteins of A thaliana and B.napus, but motifs were unevenly distributed (Figure 3). Every BnDELLA member contains four to sixteen conserved motifs, and their length ranged from six to fifty amino acids. Motif 1 to 13 were identified in all Groups except AtRGL3, BnA10RGL3, and BnC09RGL3 lacking motif 12. In which motif 7 and 8 annotated as DELLA domain (Figure S2). Moreover, Motif 14 and 15 were not detected in AtRGL2, BnA05RGL2, BnA05RGL2-2, and Group II, respectively. Motif 16 was detected in AtRGL3, BnA10RGL3, BnC09RGL3, and in BnA09GAI, BnC09GAI. Furthermore, Motif 17 was only present in the N-terminal of the AtRGA, BnC07RGA, BnA06RGA, and in BnA05RGL2, BnC05RGL2 genes. Motif 18 was detected in AtRGA and BnA09GAI, BnC09GAI, BnC07RGA, BnA06RGA. In contrast, Group I and AtRGL2, BnA05RGL2, and BnC05RGL2 had an extra motif 19 and 20, respectively. A schematic diagram of all BnDELLAs motifs logos was shown in (Figure S2). These results exhibit that the BnDELLAs subfamilies differ in motif arrangements, indicating the functional divergence of BnDELLA genes during duplication events. However, proteins with similar motifs arrangements specified the functional similarities among BnDELLA members.
Chromosome location and collinearity analysis
Chromosomal mapping determined that 10 BnDELLAs distributed on 8 B. napus scaffolds, which has not yet been assembled into a chromosome. Furthermore, no distribution of BnDELLA genes were observed in the scaffoldA01, scaffoldA03, scaffoldA04, scaffoldA07, scaffoldA08, scaffoldC01, scaffoldC03, scaffoldC04, scaffoldC05, scaffoldC06, and scaffoldC08 (Figure 4). However, six BnDELLA genes, including, BnA02RGL1, BnA05RGL2, BnA05RGL2-2, BnA06RGA, BnA09GAI, and BnA10RGL3 located on the AA subgenome. While four BnDELLAs, including BnCO2RGL1, BnCO9GAI, BnC07RGA, and BnC07RGL3 located on the CC subgenome. In contrast, BnA05RGL2 and BnA05RGL2-2 have formed a cluster in the same scaffoldA05, which might result from tandem duplication (Figure 4). These results indicate that BnDELLA genes are unevenly distributed in the B. napus genome.
Furthermore, by using the BLAST and MCScanX methods, we investigate the 6 segmental duplication pairs such as, BnA06RGA/BnC09GAI, BnC07RGA/BnA09GAI, BnA02RGL1/BnC09RGL3, and one tandem duplication BnA05RGL2/BnA05RGL2-2 were determined (Figure 5), which exhibits that during evolution segmental duplication events were the main reason for the divergence of DELLA genes in B. napus. In addition, comparative synteny of DELLA gene pairs between B. napus and A. thaliana was conducted (Figure 6). The results show that 10 BnDELLA genes of B. napus are associated with AtDELLAs, which are more than one ortholog’s copies in BnDELLAs. For instance, AtGAI and AtRGA show synteny relationships with BnA09GAI, BnC09GAI, BnA06RGA and BnC07RGA, suggesting that AtDELLAs genes might contribute to the evolution of the BnDELLAs family. Moreover, we evaluate the pressure of selective constrains on each pair of duplicated BnDELLA genes, and calculated the nonsynonymous (Ka) and synonymous (Ks) ratio (Table S2, Figure S3). Our findings showed that all of the BnDELLA pairs had the Ka/Ks ratio lower than 1, implying that BnDELLA genes experienced robust purifying selective pressure.
Prediction of the bn-microRNAs putative targets sites of BnDELLAs
The regulatory purpose of DELLAs and their interacting proteins have been characterized widely in various plant species; however, a possible underlying post-transcriptional modification of DELLAs in response to environmental stresses is still unclear [42-44]. It has been reported that miRNAs play a significant role in transcriptional and post-transcriptional level to modulate gene expression under stresses [45, 46]. To identify miRNAs interaction with BnDELLAs isoforms, we obtained the bn-miRNAs data from B. napus comprehensive study to predict the targeted BnDELLA genes sites. We found that 10 BnDELLA genes from B. napus were unambiguously complementary to the 18 known B. napus miRNAs. These miRNAs length reached from 1-24 base pairs, with 11 nt being the most frequent in all BnDELLAs (Table 2). Target prediction shows that BnDELLAsBnaC07RGA and BnaA09GAI are targeted by 2 well-known miRNAs, bna-miR6029 and bna-miR6031, respectively. Among the other bn-miRNAs identified in our analysis, bna-miR2111a, bna-miR166a are found to be involved in targeting the BnRGL1. In contrast, bna-miR172b targets BnA02RGL1 and BnA05RGL2. Additionally, bna-miR390a and bna-miR168a are found to target BnRGL3. Based on this analysis, we perceived that B. napus miRNA strongly targets B. napus both A and C genome to regulate BnDELLAs gene expression under constant stress to stabilize plant growth and defense tradeoff.
cis-element analysis in promoter regions of BnDELLAs and their distribution
Physiological and molecular studies on DELLAs suggested their role in multiple hormonal signaling pathways by interacting with a wide variety of transcriptional regulators and transcriptional factors. However, the molecular mechanism of interaction and regulation of DELLA genes are quite unclear. To gain more insights into the potential function and regulatory mechanism of the BnDELLAs, we analyzed the cis-regulatory elements in the 1500bp upstream promoter region of the BnDELLAs by using the Plant-CARE database and divided them into four categories (Figure 7A). We found that the individual DELLA gene in B. napus contains multiple cis-acting elements (Table S3), nearly all of the BnDELLA genes promoters have CAAT-box, TATA-box, light, stress, hormone, and development related responsive cis-elements. However, the distribution and numbers differed significantly between the BnDELLA genes (Figure 7B). In detail, BnA09GAI, BnA02RGL1 has a higher number of light and hormone-responsive elements. In contrast, BnA06RGA, BnA05RGL2, and BnA10RGL3 carried a higher number of stress and development related cis-elements, respectively. However, some of the cis core elements were only found in some BnDELLAs. For example, GC-motif (enhancer-like element involved in anoxic specific inducibility), DRE-core (cis-acting regulatory element regulate cold stress induce dehydration), and 3-AF binding site (part of a conserved DNA module array CMA3) were found in BnA06RGA and BnC07RGABnDELLA members. Similarly, GATA-motif (cis-acting regulatory element involved in light-responsive floral, hypocotyl, and seed development), AT-Rich sequence (cis-element for maximal elicitor-mediated activation) were present in BnCO9RGA and BnA10RGA. ATCT-motif (Part of a conserved DNA module involved in light responsiveness), Gap-Box (cis-acting element related to light-responsive GapA gene) were present in BnC02RGL1 and BnA02RGL1, respectively. Moreover, AuxRR-Core (cis-acting regulatory element involved in auxin responsiveness) was only found in BnA05RGL2 and BnA05RGL2-2. In contrast, O2-site (cis-regulatory element involved in zein metabolism regulation) was absent in all BnDELLAs except BnC09RGL3 and BnA10RGL3 (Figure 7A). These results showed that BnDELLA genes contain a wide variety of stress and defense-related cis-elements compared to development, light, and hormone-responsive cis-elements, suggesting the BnDELLAs diverse function in response to various biotic and abiotic conditions.
Transcriptomic and qRT-PCR analysis of BnDELLAs in different tissues
The candidate BnDELLA genes expression transcriptomic data from the different tissues of the B. napus variety “Zhongshuang 11” (ZS11) were obtained from the BnTIR database http://yanglab.hzau.edu.cn/BnTIR. The extracted data normalized by log2 fold change and heatmap were generated. As shown in (Figure 8), the expression pattern of the 10 BnDELLA genes was different among various tissues, which points out that the additional copies of the homologs BnDELLA genes show variations in expression during plant development. This can provide important insight into these genes’ distinct roles in B. napus. To better understand the expression pattern of the BnDELLA genes, we performed qRT-PCR in 8 primary tissues (Root, Mature-Silique, Leaf, Flower, Flower-Bud, Stem, Tips, Seed) of B. napus variant ZS11. We found a strong correlation between the transcriptomic and our qRT-PCR results (Figure 9). On the whole, BnGAI and BnRGA are highly expressed in the stem and tips, while BnRGL1 and BnRGL2 were mainly expressed in the floral organs and seed, respectively. Conversely, in our qRT-PCR analysis, BnRGL3 show minimal expression in any tissues. However, based on the transcriptomic data analysis, BnRGL3 expression was highly observed in the silique. The contradiction between the qRT-PCR and transcriptomic data, especially in the BnRGL3 expression, might be due to the harvesting of silique at 6 and 28 days after flowering, which show the complex variation of the BnDELLA genes expression at various stages of plant development. These results indicate that the expression of BnDELLA genes at multiple plant tissues might play an indispensable role in regulating gibberellins and other phytohormones signals to mediate plant growth and survival tradeoff under constant stress conditions.
Expression analysis of BnDELLAs under different stress
To further explore and gain more insights into possible BnDELLAs function. We used the publicly available RNA-seq data to analyzed the BnDELLA genes expression under different stress conditions, such as MA, CA, FA, DT, HT, ABA, NaCl, and Sclerotinia sclerotiorum. Our predicted RNA-seq data show that the BnDELLAs gene expression significantly changed upon different stress treatments. For instance, BnRGL2 was up-regulated by all denoted stresses except in drought and salt (Figure 10). Whereas, BnGAI show induced expression in response to MA, HT, DT, and salinity. In contrast, BnA10RGL3, BnC09RGL3 almost exhibits reduced expression in response to heat, drought, ABA, and salt treatment. However, induced expression was observed during cold and Sclerotinia sclerotiorum treatment. Many previous studies on AtDELLA genes had provided evidence of their distinct and fundamental role in regulating plant physiology under abiotic and biotic stress [24, 47-49], suggesting the importance of BnDELLA genes in improving stress tolerance.
In order to understand the functional regulatory mechanism, we used the AtDELLA orthologous pairs of the A. thaliana to performed GO enrichment analysis. Three common categories of GO terms were observed including, biological process (BP), cellular component (CC), and molecular function (MF). In the MF category, DELLA genes are highly enriched in binding (GO:0003700), (GO:0005515), and transcriptional regulation activity (GO:0140110). CC is enriched in the nucleus (GO:0005634), which exhibits that DELLAs are nuclear proteins. Similarly, most of the GO terms (GO:0009737, GO:0009739, GO:0009740, GO:0009753, GO:0042538, GO:0009863, GO:0072593, GO:0009651, GO:0009908, GO:2000033, GO:0030154, GO:0010187, GO:0009938, GO:0006355, GO:0010218, GO:0009723) were abundant in biological process, which show response to hormones and stresses (Figure 11, Table S7). This GO enrichment results suggested that the BnDELLAs gene plays a crucial role in regulating hormonal signaling in response to stresses, which is consistent with previous studies [9, 50-52].