Comprehensive analyses of the annexin (ANN) gene family in Brassica rapa, Brassica oleracea and Brassica napus reveals their roles in stress response

Annexins (ANN) are a multigene, evolutionarily conserved family of calcium-dependent and phospholipid-binding proteins that play important roles in plant development and stress resistance. However, a systematic comprehensive analysis of ANN genes of Brassicaceae species (Brassica rapa, Brassica oleracea, and Brassica napus) has not yet been reported. In this study, we identified 13, 12, and 26 ANN genes in B. rapa, B. oleracea, and B. napus, respectively. About half of these genes were clustered on various chromosomes. Molecular evolutionary analysis showed that the ANN genes were highly conserved in Brassicaceae species. Transcriptome analysis showed that different group ANN members exhibited varied expression patterns in different tissues and under different (abiotic stress and hormones) treatments. Meanwhile, same group members from Arabidopsis thaliana, B. rapa, B. oleracea, and B. napus demonstrated conserved expression patterns in different tissues. The weighted gene coexpression network analysis (WGCNA) showed that BnaANN genes were induced by methyl jasmonate (MeJA) treatment and played important roles in jasmonate (JA) signaling and multiple stress response in B. napus.


Results and Discussion
). All members were verified for the presence of annexin repeats using InterPro and Conserved Domain (CD)-search in NCBI. Brassicaceae species experienced an extra whole-genome triplication (WGT) event [44][45][46] , based on which approximately 24 and 48 ANN genes were expected in B. rapa/B. oleracea and B. napus genomes, respectively. However, only 13, 12, and 26 ANN genes were found in B. rapa, B. oleracea, and B. napus, respectively (Table 1). In B. napus, the number of genes in the An-subgenome (12) and Cn-subgenome (14) was almost the same as that in their diploid progenitors B. rapa and B. oleracea (Table 1). These results indicate the loss of about half of ANN genes after the Brassicaceae WGT in B. rapa and B. oleracea. However, most of the duplicated ANN genes were retained after the whole-genome duplication (WGD) event in B. napus. WGD event of gene family appears to be a widespread phenomenon, such as the auxin response factor (ARF) 47 , Auxin/Indoleacetic acid (Aux/IAA) 48 , glutathione transferases (GST) 49 , BRI1-EMS-SUPPRESSOR1 (BES1) 50 , Heat stress transcription factors (Hsfs) 51,52 , GRAS 53 family genes in diploid and allopolyploid Brassicaceae and Calcium-dependent protein kinases (CPK) 54 , Jasmonate ZIM-domain (JAZ) 55 and Nuclear factor YB (NF-YB) 56 in diploid and allopolyploid Gossypium species (G. raimondii: DD genome; G. arboretum: AA genome; G. hirsutum: AADD genome).
Among the 51 Brassica ANN genes, 35 were the typical ANN, which encoded proteins ranging from 315-325 amino acids (AA) and contained four annexin repeats. All eight ANN members (315-320 AA) homologous to AtANN4 (AT2G38750) contained 2-3 annexin repeats, as same as AtANN4. While two other ANN members (157 AA) contained only a single annexin repeat and six members (183-265 AA) contained 2-3 annexin repeats (Table 1), they may were the truncated mutated duplications.
phylogenetic and structural analysis of ANN. A phylogenetic tree (Fig. 1A) was generated using the sequences of 59 ANN proteins from B. rapa, B. oleracea, B. napus, and Arabidopsis (Fig. S1). These ANN proteins were divided into six groups. All eight AtANN were found to have orthologous genes in B. rapa, B. oleracea, and B. napus (Fig. 1A). Twelve pairs of BnaANN were found in the corresponding B. napus An-and Cn-homoeologous chromosomes, and ten pairs of them had homologous genes both in B. rapa and B. oleracea. Meanwhile, two pairs (BnaC03g49290D/BnaA06g23960D and BnaC09g44350D/BnaA10g20320D) only had homologous genes in B. rapa. All 12 BoANN genes were found to have homologous genes in the Cn-subgenome of B. napus, while one BrANN (Bra039578) had no homologous gene in An-subgenome of B. napus (Table 1 and Fig. 1A).
Gene structure analysis revealed that majority of the homologous ANN gene pairs had same gene structure (Fig. 1B). There were five introns in group IV/V/VI members, expect for two truncated mutant genes (Bo1g039570 and BnaC01g16910D) (Fig. 1B). This finding indicates that the ANN genes are conserved in Brassicaceae species, possibly due to their importance in plant growth and productivity.
A typical ANN protein contains four annexin repeats, each approximately 70 amino acids long 1,3 . Annexin repeat usually contain a characteristic type II motif for binding calcium ions with the sequence GxGT-[38 residues]-D/E 3 . MEME analysis showed that 42 ANN proteins contained four annexin repeats (Fig. 1C). Motif1 was the core sequence of all the four annexin repeats, and motif4 was only found in the third annexin repeat in group I-V, while motif5 was the core sequence close to the C-terminal of Motif1 in the second and fourth annexin repeats (Fig. 1C). According to the gene structure and motif analysis, the missing parts of the truncated mutant members were readily apparent. Both the first and fourth annexin repeats were absent in Bo9g172330 and BnaC09g46400D, and the first annexin repeat was absent in BnaAnng37420D. Bo1g039570 and BnaC01g16910D had only the second annexin repeat at the C-terminal (80-159 AA), and the core sequence of annexin repeat was not detected at the N-terminal (1-79 AA). It is similar in the N-terminal of Bo6g043900, BnaC06g08410D, and AtANN4 homologues ( Fig. 1B,C). chromosomal location and synteny analysis of ANN of B. rapa, B. oleracea, and B. napus. As showed in Fig. 2, the distribution of BnaANN in An-and Cn-subgenome was nearly even with 12 ANN genes from the An-subgenome and 14 from the Cn-subgenome. However, the ANN genes' distribution was uneven on each chromosome. Three pair (2 genes/pair) of ANN genes from the An-subgenome were repeated in tandem on chromosome Bn_A03, Bn_A04, and Bn_A10 (Fig. 2); and three pair (2 genes/pair) of ANN genes from the Cn-subgenome were repeated in tandem on chromosome Bn_C03, Bn_C04, and Bn_C09 (Fig. 2C). B. napus genome analysis showed that the An-and Cn-subgenome were largely collinear to the corresponding diploid Ar and Co genomes 43,57 . Most of the An-Ar and Cn-Co orthologous gene pairs demonstrated similar chromosomal locations. The distribution of ANN genes in B. rapa and B. oleracea were similar to the distribution of the orthologous BnaANN genes in the B. napus An-subgenome and Cn-subgenome, respectively (Fig. 2). Two BnaANN (BnaAnng04520D and BnaAnng37420D) and one BrANN (Bra039578) genes were located on the unanchored scaffolds that were not mapped to a specific chromosome (Fig. 2). The sequence and phylogenetic analyses revealed BnaAnng04520D-Bra034402 and BnaAnng37420D-Bra031890 as two An-Ar orthologous gene pairs. Based on this, we predicate Bn_A02 and Bn_A05 as the true chromosomal locations of BnaAnng04520D and BnaAnng37420D, respectively. BnaANN (BnaC03g49290D and BnaC09g44350D) had no orthologous genes in B. oleracea (Fig. 2), though they had homologous genes in An-subgenome. These findings indicate that duplication of BnaA06g23960D and BnaA10g20320D led to the formation of BnaC03g49290D and BnaC09g44350D, respectively. Analysis of the synteny among An-subgenome and Cn-subgenome showed high collinearity between Bn_ A01-Bn_C01, A02-C02, A03-C03, A04-C04, A05-C05, A06-C06, A07-C07, A08-C08, A09-C09, and A10-C09, and 83.7% orthologous gene pairs between B. rapa and B. oleracea were retained as homologous gene pairs in B. napus 43,57 . 90.9% ANN gene pairs (10/11 pairs) between B. rapa and B. oleracea were retained as homologous gene pairs between B. napus An-chromosomes and Cn-chromosomes (Fig. 2). www.nature.com/scientificreports www.nature.com/scientificreports/ There were two tandem pairs (AtANN3/4 and AtANN6/7) on chromosome 2 and chromosome 5 in Arabidopsis, respectively 58 . Bra009048/Bra009049, Bo9g172330/Bo9g172340, BnaA10g22010D/BnaA10g22020D, and BnaC09g46400D/BnaC09g46410D were homologous to AtANN6/7 tandem pair in B. rapa, B. olearcea, B. napus An-subgenome and Cn-subgenome, respectively. We identified two tandem pairs each homologous to AtANN3/4 in B. rapa (Br_A03 and Br_A04), B. olearcea (Bo_C03 and Bo_C04), B. napus An-subgenome (Bn_ A03 and Bn_A04), and Cn-subgenome (Bn_C03 and Bn_C04) (Fig. 2). AtANN8 (AT5G12380) was located near the AtANN6/7 tandem pair on chromosome 5 in Arabidopsis 58 . Correspondingly, there was a gene homologous to AtANN8 located near the tandem pair homologous to AtANN6/7 in B. rapa (Br_A10), B. napus An-subgenome (Bn_A10) and Cn-subgenome (Bn_C09) (Fig. 2). There was no gene homologous to AtANN8 in B. olearcea (Bo_C09). Instead, we found a truncated mutated gene (Bo1g039570) homologous to AtANN8 in B. olearcea (Bo_C01). Meanwhile, a truncated mutated gene (BnaC01g016910D) was homologous to Bo1g039570 in B. napus Cn-subgenome (Bn_C01) (Fig. 2).   (Table S2) 57,60,61 . The ANN genes were expressed across different vegetative and reproductive organs during different developmental stages of the four species (Fig. 3). In general, the ANN expression pattern was different between groups; however, expression pattern was very similar within a group in the four plant species.
Group I (ANN6/7) members showed expression in young siliques (ovules) and seeds, which indicate their importance in ovule and seed development in Brassicaceae plants. However, two truncated mutated members (Bo9g172330 and BnaC09g46400D) homologous to ANN7 were at low abundance expression levels (Fig. 3). Unlike Bo9g172330 and BnaC09g46400D, other five truncated mutated members (BnaAnng37420D, BnaA06g23960, BnaC06g08410D, Bo1g039570 and BnaC01g16910D) have a similar expression level to their homologous genes which have complete gene structure (Figs. 1 and 3). So, truncated mutated gene structures may decrease their own genes' expression level, but not always. The expression levels of group 2 (ANN2) members were highest in roots and young siliques (ovules), while that of group 3 (ANN1) members were higher in roots, stems, and young siliques (pericarps) (Fig. 3). These expression levels are consistent with the role of AtANN1 and AtANN2 in root growth and development [20][21][22] . It was indicated that ANN1/2 regulates the development of young siliques and seeds. We detected low level of expression for group 4 (ANN8) members. AtANN5, which regulates pollen development 23,24 , showed specific expression in mature pollen. The B. napus genes homologous to AtANN5 were mainly expressed in buds and new pistils. The genes homologous to AtANN3 and AtANN4 demonstrated similar expression pattern. Both genes were expressed in flowers and young siliques (ovules), though they belong to group IV and VI, respectively (Fig. 3). All these indicated that ANN genes may be involved in various developmental processes with different functions. In Arabidopsis, AtANN3 and AtANN4 had similar expression pattern because they share a 5′ promoter region (2654 bp) 58 . In B. rapa, Bra000090 and Bra000091 share a 5′ promoter region (2079 bp), while in B. oleareca, Bo3g032760 and Bo3g032770 share a 5′ promoter region (2412 bp); In B. napus, BnaA03g18070D and BnaA03g18080D share a 5′ promoter region (2455 bp) and BnaC03g21590D and BnaC03g21600D share a 5′ promoter region (6199 bp). They were homologous to AtANN3/AtANN4 pair, and had similar expression pattern. But another gene pairs (Bra017102/Bra017103, Bo4g187790/Bo4g187780, BnaA04g22190D/BnaA04g22180D, and BnaC04g45920D/BnaC04g45910D) homologous to AtANN3/AtANN4 pair were at low abundance expression levels (Fig. 3B,C). All the results suggested that there were gene duplications, gene expression pattern differentiations and subsequent functional diversifications in ANN family genes in Brassicaceae species, and the functions of homologs of a given group ANN genes might be redundant as they share similar expression patterns. Accumulating evidence from various plant species has shown the regulation of ANN genes in response to abiotic stress and hormonal treatment [5][6][7]9,58 . To examine the expression pattern of BnaANN genes under various abiotic stress conditions and hormonal treatments, we utilized the data on transcriptional profiling (Table S3). As shown in Fig. 4, most of the expressed BnaANN genes in group II/III/V/VI/were up-regulated under salinity and PEG stress in roots and MeJA treatment in leaves. BnaA06g23960, BnaA03g18070D and BnaC03g21590D were down-regulated under cold stress, whereas BnaAnng04520D and BnaC05g27530D were up-regulated under cold stress at 12 hours point (Fig. 4).
B. napus is a winter biennial crop with excellent tolerance to low-temperature stress during vegetative stage. The response mechanisms are different under chilling and freezing temperatures, as well as cold shock and cold acclimation in plants 62,63 . Based on the transcriptional profiling of early-maturing, cultivated B. napus varieties under different low-temperature treatments with or without cold acclimation (GSE129220: https://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE129220) (Table S4) 64 , transcriptome analysis revealed that group III BnaANN were induced slightly by chilling stress, and were up-regulated by freezing stress strongly, regardless of cold acclimation (Fig. 5A). This finding indicates that group III BnaANN genes play important roles in freezing stress in B. napus.
Sclerotinia sclerotiorum is a hemibiotroph pathogen with a wide host range. It is the causative agent of stem rot, one of the most devastating diseases of B. napus 65,66 . Previous studies have shown the role of JA signaling in plant resistance to hemibiotroph pathogens [67][68][69][70] . The transcriptional profiling of B. napus susceptible (Westar) and tolerant (ZY821) genotypes infected with S. sclerotiorum (GSE81545: https://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE81545) (Table S5) showed that the group II BnaANN were induced by S. sclerotiorum infection, and the expression level in the susceptible genotype (Westar) was more than that in the tolerant (ZY821) genotype; some members from group III and group V BnaANN were induced, while some members were repressed by S. sclerotiorum infection (Fig. 5B). These findings indicate a complex response mechanism and the role of some BnaANN in B. napus response to S. sclerotiorum.
To validate the results on transcriptional profiling, we performed a qRT-PCR to detect the transcript levels of three genes (BnaC03g49290D, BnaC05g27530D, and BnaC03g21590D) from group II/III/VI in the roots challenged with salt and PEG and in the leaves treated with cold and MeJA. The expression pattern (of these BnaANN genes) was consistent with the RNA-Seq data (Figs. 4-6). All three BnaANN genes were induced under salinity and PEG stress in roots and induced by MeJA in leaves (Fig. 6A). BnaC03g49290D and BnaC03g21590D were repressed by cold treatment (Fig. 6A). BnaC05g27530D was significantly upregulated under freezing stress, with or without cold acclimation (Fig. 6B). In B. rapa, Bra034402 (gene to homologous BnaC05g27530D) was strongly induced by hormone and stress treatments 11 . All these results indicated the role of these three genes in multiple abiotic stress response and JA signaling response in B. napus.  Table S6). The lightgreen module (845 genes) was positively correlated with the MeJA treatment in leaves (Fig. S2). Five  BnaANN genes (BnaA03g18070D, BnaA03g18080D, BnaC03g21590D, BnaC03g21600D and BnaC05g27530D) were induced by MeJA treatment in lightgreen module ( Fig. 4 and Table S6). The top two hub genes with the highest the module membership kME (k-means clustering algorithm) values were BnaA03g18070D (BnaANN4A-1) and BnaC06g31830D (BnaTIFY7) in the light green module (Fig. 7A and Table S6). The jasmonate acid (JA) signaling repressor, TIFY, was induced by JA and regulates plant development and stress response [72][73][74][75] . Additionally, there were some B. napus JA biosynthesis genes and JA responsive genes in the light green module, such as the Lipoxygenase (LOX), Allene oxide cyclase (AOC), Allene oxide synthase (AOS), 12-oxophytodienoate reductase (OPR), Jasmonate O-methyltransferase (JMT), and Ethylene-responsive factor (ERF) (Fig. 7A and Table S6). Transcriptional profiling and qRT-PCR analysis results showed that BnaA03g18070D/BnaANN4A-1, BnaC06g31830D/BnaTIFY7, BnaC04g38070D/BnaERF42, and BnaC02g29610D/BnaAOS were all induced by MeJA ( Fig. 7B and Table S7). However, there was little research at the functions of annexins in JA signaling. ZmAnx6.1 and ZmAnx7 were induced at 12 h by JA, and the JA-responsive cis-elements exist in their promoters 76 . We analyzed the promotor sequences (2000 bp upstream of transcription start sites) of BnaANN, and founded that there were so many cis-elements involved in stress (drought, low-temperature, heat, anaerobic, wounding, defense and stress) response and plant hormones (MeJA, ABA and SA) response in their promotors, MeJA-responsive cis-element (CGTCA-motif, TGACG-motif and G-box) was the most numerous cis-element and all BnaANN members contain MeJA-responsive cis-elements (1 to 9) in promotors (Fig. S3). It suggested that the BnaANN genes in lightgreen module involved in JA signaling response in B. napus.
Three BnaANN genes (BnaA04g22190D, BnaC03g49290D, and BnaC02g43450D) in the blue module were expressed with NaCl and PEG treatments in roots, while genes (BnaC08g06690D, BnaA10g20320D www.nature.com/scientificreports www.nature.com/scientificreports/ and BnaC09g44350D) in the green module were expressed in roots. BnaC01g16910D, BnaA06g23960D, and BnaC06g08410D in the turquoise module were positively correlated with bud, stamen, ovule, and silique ( Fig. S2 and Table S6). All the results indicate the different functions of B. napus ANN genes during plant development and stress response. ANN of B. rapa, B. oleracea, and B. napus. B. rape, B. oleracea and B. napus ANN proteins have been identified using BLASTP (E-value < 1e-5) to look for homologs of Arabidopsis ANN among B. rape, B. oleracea and B. napus genome sequences database in Ensembl gemones (http://ensemblgenomes.org/) 77 . The annexin motifs in ANN proteins were characterized using InterPro (http://www.ebi.ac.uk/interpro/) 78 and the NCBI conserved domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi).

Identification of
The molecular weight (Mw), isoelectric point (pI), and subcellular localization of ANN proteins were predicted using the Compute pI/Mw tool (http://web.expasy.org/compute_pi/) 79 and ProtComp 9.0 (http://linux1. softberry.com/). The exon and intron organization of the ANN genes were analyzed using the Gene Structure Display Server (GSDS) (http://gsds.cbi.pku.edu.cn/) 80 . The conserved motifs of ANN were analyzed with MEME (http://meme.nbcr.net/meme/cgi-bin/meme.cgi) 81 . nonsynonymous and sunonymous substitution rate ratio (Ka/Ks). DnaSP (DNA Sequence Polymorphism) v6 84 was used to calculate the ratio of the nonsynonymous substitution rate (Ka) to the synonymous substitution rate (Ks) and the Ka/Ks value between paralogous gene pairs. plant materials and treatments. ZS11 (B. napus L. cv. Zhongshuang 11) 57 seeds were allowed to germinate and then the seedlings were transplanted to pots containing soil or vermiculite. The growth conditions, hormone treatments, and abiotic stress conditions were as described previously 85 . Hormone treatments were performed by spraying leaves of 6-week-old seedlings with ABA (100 μM), MeJA (100 μM), SA (1 mM), and ETH (10 μg/ml); To simulate hot and cold stresses, seedlings were grown in chamber with 40 °C or 4 °C. To simulate salt and PEG stresses, seedlings were irrigated with NaCl (200 mM) or PEG-6000 (20%) solutions.
For chilling and freezing treatments with or without cold acclimation, the seedlings of two early-maturing semi-winter rapeseed varieties (HX17 and HX58) were used. They were treated as described previously 64  RnA isolation and sequencing and gene expression analysis. The collected samples were sent to the sequencing cooperations of Sangon Biotech (Shanghai) Co., Ltd. and Novogene Co., Ltd. for RNA isolation, examination, and sequencing 64,85 . qRT-PCR analysis was performed as described previously 85 . The primers used in this study were listed in Table S8.