Systematic analysis of JmjC gene family of Brassica napus and KDM5 subfamily involved in abiotic stress response

Background: Jumonji C (JmjC) proteins play an important role in plant development and stress response through the removal of lysine methylation from histones. Brassica napus, which originated from spontaneous hybridization between Brassica rapa and Brassica oleracea, is the most important oilseed crop after soybean, but evolutionary relationships and functions of JmjC proteins remain unclear. Results: 65 JmjC genes were identified from B. napus genome, 29 from B. rapa, and 23 from B. oleracea. These genes were grouped into seven clades according to conserved sequences, and their catalytic activities of demethylation were predicted. Group-KDM4/JHDM3 for H3K4/9/27/36, Group-KDM5A/B for H3K4, Group-JmjC domain-only A/B for H3K27/36, Group-KDM3/JHDM2 for H3K9, and Group-JMJD6 may be for arginine demethylases. B. napus inherited most of its JmjC genes from its parents. The average retention rate of B. napus JmjC gene from B. rapa (93.1%) and B. oleracea (82.6%) exceeded that of all homologous gene pairs (83.7%) across the whole B. napus genome. Thirteen new or duplicated JmjC genes have emerged in B. napus. Sequence similarity and domain organization analyses suggest that the functions of these genes might be diversified. Furthermore, KDM5 genes were examined under stress conditions due to H3K4 demethylation. Expression profiles indicated that the genes from B. napus are possibly involved in various stress responses. Conclusion: This study provides the first genome-wide characterization of JmjC genes in Brassica species. Its JmjC genes potentially have diverse functions, and its KDM5 genes might be involved in stress response. The results of this study facilitate the future functional characterization of the demethylation of in

3 residues, form "histone code" to regulate various biological processes [2]. Histone methylation is usually catalyzed by three protein families of histone methyltransferases: protein arginine methyltransferase family, Su (var)3-9, Enhancer-of-zeste and Trithorax (SET) domain family, and telomeric silencing disruptor that is also known as DOT1-Like (Kmt4/DOT1L) [3]. Histone lysine methylation, playing many different roles in biological processes ranging from heterochromatin formation to transcription regulation, is dynamically regulated by histone lysine methyltransferases (KMTs) and histone lysine demethylases (KDMs), and can be distinguished depending on the position of lysine residue and the number of added methyl groups in lysine residues, which carry mono-, di-, or tri-methylated groups [4].
As sessile organisms, plants are more susceptible to environmental factor than animals, though the former has developed complex mechanisms that perceive environmental stimuli and respond and adapt accordingly to the changing external environment. Histone modifications play important roles in a wide range of biological processes by affecting gene expression. For example, H3K4me3 and H3K36me3 are correlated to active gene expression, and H3K9me2 and H3K27me3 are associated to gene silencing [5]. The reversible epigenetic modifications dynamically regulate the plant's response to stress, including cold, freezing, saline, drought and submergence [6]. Abiotic stresses might induce various histone modifications to regulate the expression of genes with enhancing stress tolerance [2].
The active change in histone lysine methylation denotes a crucial role in plant adaptation to adverse environmental conditions. Athaliana genome-wide H3K4 methylation patterns (H3K4me1, H3K4me2 and H3K4me3) show dynamic responses to dehydration stress [7]. Arabidopsis ATX1, H3K4me3, is involved in dehydration response through ABA-dependent and -independent pathways [8]. Under drought stress, H3K4me3 enrichment is correlated with the activation of Arabidopsis drought stressresponsive genes, such as RD29A and RD20 [9,10]. H3K4me3 can be maintained at low levels after rehydration and which could function as an epigenetic mark of drought stress memory [10]. Under heat stress, H3K4 methylation accumulates to induce gene expression and can be sustained after heat stress to positively respond on a future stress incident [11]. High-dose oxidative stress increases several histone methylation marks, including H3K4me3, and low-dose oxidative stress increases the global levels of H3K4me3 and H3K27me3 [12]. Rice genome-wide H3K4me3 profiling showed that H3K4me3 modification is remarkably and positively correlated with the transcript level of drought stress-responsive genes [13]. Under submergence stress, the H3K4 methylation status of rice submergence-inducible ADH1 and PDC1 genes switches from dimethylation to trimethylation at the coding regions and is correlated with the increased expression of ADH1 and PDC1 [9].
Histone methylation is also a reversible process regulated by methylases and demethylases and might involve three distinct classes of demethylase. Petidylarginine deiminase 4 was identified as the first enzyme capable of antagonizing histone arginine methylation, but it cannot be strictly considered as a histone demethylase for its target citrulline [3]. KDMs mainly consist of LSD1/KDM1s (Lysine specific demethylase 1) and JmjC-domain enzymes, which both utilize oxidative mechanisms. LSD1, a flavin-dependent lysine specific demethylase member, was first reported as a histone lysine demethylase and transcriptional co-repressor specifically for histone H3 lysine 4 in humans [3].
Arabidopsis homologues of human LSD1 also exhibit similar demethylation and reduce H3K4 methylation in chromatin to promote floral transition [14].
Allotetraploid species Brassica napus (oilseed rape, AACC, 2n=38) originated from interspecific spontaneous hybridization between Brassica rapa(AA, 2n = 20) and Brassica oleracea(CC, 2n=18) [21]. Brassica species might have diverged from a common ancestor with an Arabidopsis lineage from 14.5-20.4 million years ago [25]. The protein organization and function of JmjC domain in Brassica species and its relative relationship with model plant Arabidopsis remain uncharacterized. B. napus is currently the most important oilseed crop, preceded only by soybean. However, B. napus is vulnerable to abiotic stress that limits its growth and productivity and reduces its economic benefits.
KDM5/JARID1 subfamily may regulate many abiotic stress responses genes through H3K4me3 and H3K4me2 downregulation but the roles of H3K4 demethylation in abiotic stress remain unknown.
The JmjC genes in A and C subgenomes of B. napus show nearly identical distributions to its ancestor genomes B. rapa (A-genome, 29) and B. oleracea (C-genome, 23) (Fig. 1). A02 and A07 chromosomes only exist in one member of B. napus, which is similar to its ancestor B. rapa genomes. A09 chromosome carries the highest number, seven genes. Four tandem JmjC genes pairs located on chromosomes A03, A09, and C03 in B. rapa (Fig. 1)  domain is the second highly-conserved domain that is close to the N terminus and shorter than JmjC domain [16]. The four tandem array ZnF domain of RELATIVE OF EARLY FLOWERING 6 (REF6)/AtJMJ12 targets motif CTCTGYTY, and the ZnF domain only exists in subgroup-I [27]. REF6 also tends to bind to hypo-methylated CTCTGYTY motifs in vivo [28]. Subgroup-I generally harbors 7-8 exons, but subgroup-II keeps highly similar gene structures with 10 exons (Fig. 3C).
JmjC proteins have been discovered as Fe(II)-and αKG-dependent histone demethylases [17,18]. The JmjN and JmjC domains, two non-adjacent domains, interact with each other through two ß-sheets and form a single functional unit to ensure the stability and appropriate transcription activity of Gis1 and maintain the overall protein levels and function of Jhd2 H3K4-specific demethylase in budding yeast [29,30]. KDM4/JHDM3 has conserved Fe(II) binding site (His and Gluresidues) and αKG binding site Group-KDM5B differs from group-KDM5A group in domain organization, which has BRIGHT and PHD but lacking FYRN and FYRC (Fig. 4F). BRIGHT is associated with H3K4 demethylase by DNA binding motif (CCGCCC) to regulate transcription [33]. PHD mainly exerts epigenetic effectors capable of recognizing or "reading" post-translational histone modifications and unmodified histone tails [34].
The original PHD role in gene transcription is acted as a reader of H3K4me3 in 2006 [35]. Many sophisticated functions of PHD were also determined, including H3K9me3 recognition and binding to the N-terminus of H3, indicating its key roles in regulating transcription and chromatin structure [36]. Group-KDM5A/B shows a wide range intron/exon number , but sister gene pairs are relatively conserved in gene structure (Fig. 4). In group-KDM5A, subgroups-I/Ⅱ are highly conserved in Fe(II) and αKG binding sites, except for BrJMJ16;b in which Phe is replaced by Met in αKG binding site, and His is replaced by Arg Fe(II) binding site. In subgroup-III, Phe is replaced by Gln in αKG binding site, and BoJMJ19;c/d is variable in other Fe(II) and αKG binding sites (Fig. 4D). In group-KDM5B, BnJMJ17a gene structure is similar to its parent BrJMJ17 (Fig. 4G). Group-KDM5B is highly conserved in Fe(II) and αKG binding sites, similar to KDM4/JHDM3 group (Fig. 4H). Group-JmjC domain-only A has stable exon distribution harboring approximately 7-9 exons. Sequence alignment and logos analysis of JmjC domain reveal that JmjC domain-only A group is highly conserved in Fe(II) and αKG binding sites. However, compared with that in the KDM4/JHDM3 group, Phe is replaced by Thr in Fe(II) binding site (Fig. 5D). In group-JmjC domain-only B, subgroup-I genes contains 6 exons, subgroup-III harbors 4 exons, and subgroup-II has many exons (Fig. 5G). As compared with that in the KDM4/JHDM3 group, the Phe residue is replaced by Ser within AtJMJ31 orthology (Fig. 5H).

Group-JMJD6
The phylogenetic tree showed that the JMJD6 group is close to JmjC domain-only A group and includes  (Fig. 6). F-box domain recognizes a wide array of substrates and regulates many important biological processes by degrading cellular proteins in plants [37].
In group-KDM3/JHDM2 (Fig. 7B), half of the member harbors RING domain as the second primary domain. Cys-X 2 -Cys-X (9-39) -Cys-X (1-3) -His-X (2-3) -Cys-X2-Cys-X (4-48) -Cys-X 2 -Cys is the canonical RING [38]. The RING domain of many proteins mainly binds to ubiquitin-conjugating enzymes and mediates the direct transfer of ubiquitin to substrate [38]. The AT-hook is a small DNA-binding motif with a preference for A/T rich regions found in various proteins, such as the high mobility group proteins [39].  (Fig. 2). The average retention rates from ancestor exceed the rate of all homologous gene pairs (83.7%) across the whole B. napus genome [43]. Each member of B. rapa and B. oleracea can be paired to at least one 13 homologue of B. napus, except for five members of KDM5A subfamily: BrJMJ14;b, BrJMJ15;b, BoJMJ15;b, BoJMJ19;c and BoJMJ19;d. This finding indicated their absence during allotetraploid formation (Fig. 4). Gene duplication expands genome content and changes gene function to ensure the optimal adaptability and evolution of plants [44]. The 65 JmjC proteins from B. napus were more than the total number of proteins for B. rapa (29) and B. oleracea (23) (Addition file 1). According to the systematic analysis of JmjC proteins (Figs. 2-7 These genes pairs were duplicated through segmental duplication (Addition file2). The parent of BnJMJ17;c was not found using the method by Yang et al. [45] and Sun et al. [46], but the JmjC domain sequence was consistent with the BnJMJ17;b.

Conservation and function of JmjC proteins of B. napus
Histone modification regulates plant development events by epigenetically silencing or activating target gene expression. Histone methylation is an important method to control plant development and stress response. In general, H3K4me2/me3 and H3K36 correlates with transcriptional activation, and H3K9me2 and H3K27me3 correlates with gene silencing [5]. Histone modification status is regulated by KMTs and KDMs. KDMs can balance histone methylation through the removal of methylated residue of histone lysine or arginine. Histone demethylation, catalyzed by JmjC proteins, mainly occurs at the Lys residues of histone H3 (K4, K9, K27, and K36), which is involved in plant developmental stages [47], defense against [48,49], proteasomal degradation [50] and circadian clock [51]. In addition, JmjC protein regulates cell cycle progression by regulating H4K20 methylation [52]. Recent studies showed that arginine demethylation can also be catalyzed using a subset of JmjC KDMs [53].
Arabidopsis H3K4me3 of AHG3, catalyzed by ATX4 and ATX5, plays an essential role in drought stress response [73]. Arabidopsis H3K4 hypermethylation is associated with transcriptional activation and maintenance in heat stress response [74].
The expression patterns of BnKDM5 show its involvement in drought, high temperature and salt stress 16 response (Fig. 8). KDM5B subgroup have three copies of AtJMJ17 (BnJMJ17;a/b/c). The similar gene structures and domain organizations of JmjC genes between Arabidopsis and B. napus suggest their conserved biological functions (Fig. 4). Under drought, high temperature or salt stress, all members of BnKDM5B exhibited remarkable elevated expression, except for BnJMJ17;a under salt stress, suggesting that these homologous genes have similar responses to similar stress stimuli. Moreover, some members of BnKDM5B were relatively conserved between their homologous genes in certain stress. For instance, the homologous gene of AtJMJ15 (BnJMJ15;a/b/c), AtJMJ16 (BnJMJ16;a/b/c/d), and

Domain organization and phylogenetic analysis
The gene structures were visualized using the Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/). The site information of the domain organization was used to construct a protein organization sketch map using DOG program [76].
Multiple sequence alignment of JmjC proteins ClustalW [77] and its resulting file were subjected to phylogenic analysis using MEGA7.0 program [78]. A tree was constructed based on the full-length protein sequences using neighbor-joining method with pairwise deletion and p-distance model, and a Bootstrap test of 1000 replicates for internal branch reliability. CDS sequences of JmjC were aligned with ClustalW (http://www.genome.jp/tools/clustalw/), and the resulting files were used to create Logo maps (http://weblogo.berkeley.edu/logo.cgi).

Duplicated JmjCgenes in B. napus
Duplicated genes were defined according to Yang et al [45] and Sun et al [46]. The full-length CDS sequence coverage and amino acid identities were determined using Blastn/Blastp at the NCBI website. The number of non-synonymous mutations (Ka) and the number of synonymous substitutions (Ks) of duplicated genes were calculated by DnaSP 6.0 [79]. The Ka/Ks ratios between duplicated genes and gene pairs were analyzed to determine the mode of selection. The duplication time (T, million years ago, MYA) was calculated as T = Ks/2λ× 10 −6 (λ = 1.5 × 10 −8 ) for B. napus [80]. For salt stress, seedlings were treated with 0, 100, 200, and 300 mM NaCl, and leaves were harvested 3 days after treatment. For high temperature stress, seedlings were grown in 40 °C, and leaves were harvested at 0, 12, 24, and 36 h after treatment. All harvested samples were immediately frozen in liquid nitrogen. Three independent biological replicates for each treatment were conducted.
RNA extraction and real-time quantitative PCR (RT-qPCR) RNA extraction and real-time quantitative RT-PCR were conducted as described previously [81]. Total RNA was extracted from samples using a TRIzol reagent kit (Invitrogen, Carlsbad, CA, US) following manufacturer's instructions. Total RNA were reverse transcribed into cDNA using Revert Aid RT Kit The field trail experiments in the current study were permitted by the local government in China.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Availability of data and material
All data generated or analyzed during this study are included in this published article and its Additional files. The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.   using ClustalW, and the phylogenetic tree analysis was performed using MEGA7.0.The trees were constructed with the following settings: tree inference as neighbor-joining; include sites as pairwise deletion option for total sequences analysis; substitution model as p-distance.         constructed with the following settings: tree inference as neighbor-joining; include sites as pairwise deletion option for total sequences analysis; substitution model as p-distance.