CBF transcription factors involved in the cold response of Camellia japonica (Naidong)

Identification of new CjCBFs genes in C. japonica (Naidong) involved in cold tolerance

We used our C. japonica (Naidong) low temperature and drought full-length transcriptome data based on PacBio and Illumina sequencing (PRJNA689105) as the local database. The protein sequences of CBF family genes (i.e., AtCBF1, AtCBF2, AtCBF3, AtCBF4) in Arabidopsis thaliana were queried to search the local databases using the BLASTP method with an E-value threshold of <1e−20. Then, HMMER3.0 (Potter et al., 2018) was used to further screen sequences using default parameters (E<1e−5), the AP2 domain was based on Pfam database (Pfam accession, PF00847). The longest protein was selected as a representative for the gene. Sequences that do not contain the complete AP2 domain and signature sequences were eliminated.

Sequence analysis

ExPASy (Gasteiger et al., 2003) was used to analyze the relative molecular weight, isoelectric point, hydrophilicity, and other basic information of CjCBFs according to default parameters.

Motif and phylogenetic analysis

MEME (http://meme-suite.org/tools/meme) was used for the conservedmotif analysis of CjCBFs and AtCBFs proteins, with the maximum number of motifs set at ten. The parameters were as follows: minimum motif width, 6; maximum motif width, 50; Selected the motif classic mode and 0-order model of sequences. MAFFT was used for sequence alignment according to default parameters (FASTA format, auto). The R package ggmsa was used to visualize the sequence alignment. In this study, C. japonica (Naidong) and 18 other plant species were selected from dicotyledons and monocotyledons, including Brachypodium distachyon, Oryza sativa and three additional monocotyledons, as well as Camellia japonica, Arabidopsis thaliana, Cucumis sativus and 19 other dicotyledons. The CBF-like protein data of the 19 species (with the exception of C. japonica) were downloaded from the NCBI Database. The protein sequences were selected by HMMER3.0. If the AP2 domain was truncated, or the AP2 domain match E-value exceeded 1e−5, the protein sequences were excluded. Based on this, 424 CBF-like genes were identified from the 19 plant species in order to construct the phylogeny tree. The phylogeny was constructed using the maximum likelihood method with 1,000 bootstrapping replicates in MEGA (v. 6.0) (Tamura et al., 2013). The R package ggtree (v2.0.4) was used to visualize the phylogenetic tree.

GO and interaction analysis

KOBAS (Xie et al., 2010) (http://kobas.cbi.pku.edu.cn/) and A. thaliana (thale cress) were used for the GO enrichment analysis. The sequences of four CjCBFs were inserted separately into pGADT7 Vector and pGBKT7 Vector; pGADT7-T+pGBKT7-Lam was the negative control and pGBKT7-53+pGADT7-T was the positive control. In order to verify self-activation, the pGADT7+ pGBKT7 - CjCBF1/ CjCBF2/ CjCBF3/ CjCBF5 plasmid was transformed into Y2HGold yeast competent cells, with 30 mg/ml α-gal employed to detect the self-activation. Furthermore, in order to determine the suitable concentration of AbA, the inhibitory self-activated AbA concentrations of 0.20 µl/ml, 0.30 µl/ml, 0.50 µl/ml, 0.70 µl/ml and 0.90 µl/ml were adopted for the five working solutions. The pGBKT7 and pGADT7 + CjCBF1, pGBKT7 and pGADT7 + CjCBF2, pGBKT7 and pGADT7 + CjCBF3, and pGBKT7 and pGADT7 + CjCBF5 plasmids were transformed into Y2HGold yeast competent cells, respectively, to exclude false positives. Moreover, the pGADT7+ CjCBFs and pGBKT7 + CjCBFs plasmid was transformed into Y2HGold yeast competent cells to verify the interaction between CjCBFs.

RNA-seq data

In order to explore the transcriptional expression of CjCBFs at low temperature stress across time, we evaluated our low temperature stress transcriptome (PRJNA689105) under a 4 °C treatment lasting 0 h, 12 h, 24 h, and 72 h, with four groups of transcriptome data. The expression of the CjCBFs gene was analyzed and the results were visualized by the R package pheatmap. qRT-PCR technology was employed to verify the transcriptome data and the co-expression method was used to predict the potential lncRNA regulating CjCBFs (Li et al., 2015).

Subcellular localization

Four DNA sequences encoding CBF proteins were inserted into the super1300-GFP overexpression vector, and Agrobacterium tumefaciens GV3101 containing super1300-CjCBF1, super1300-CjCBF2, super1300-CjCBF3, and super1300-CjCBF5 was injected into the Abaxial of tobacco (Nicotiana benthamiana) leaves (Li et al., 2013). A Leica DM2500-DM2500 LED fluorescence microscope was used to obtain fluorescence information.

Promoter cloning and Cis-Element analysis

C. japonica (Naidong) DNA was extracted via the CTAB method (Doyle & Doyle, 1987) to clone the promoter. The CjCBF1 upstream 5′ UTR sequence was obtained from the Takara Genome Walking Kit (Takara, Japan) following the manufacturer’s instructions. Design of nested PCR specific primers with primer 5.0 (Table S7). The obtained sequence was connected to a T vector for sanger sequencing and the Cis-Elements in the sequence were analyzed based on the PlantCARE database (Lescot et al., 2002).

Plant materials and low temperature treatment

C. japonica (Naidong) seeds were obtained from the Qingdao Botanical Garden (36°05′N, 120°08′E), Qingdao Agricultural University Cooperative, China, in September 2018. The seeds were sown in sand in the Qingdao Agricultural University culture room (20 ± 2 °C) during the winter until the fourth true leaf appeared. Following germination, the plants were transferred into the artificial climate room of the university (36°31′N, 120°39′E) under normal culturing conditions (plastic pots: top/bottom diameter 20/12 cm, height 26 cm; temperature: 25 ± 2 °C; relative humidity: 60%; soil: peat soil 80% + river sand 20%; natural light) for 3 months.

For the tissue specific expression experiment, we selected seedlings of a similar growth status and divided them into two groups, with three plants per group. Each group was treated at 4 °C for either 0 h or 24 h in a low temperature incubator. Tender roots, stems, leaves, buds, and flowers were collected at 0 h and 24 h. For the heat, cold and light treatments, seedlings were separated into several groups. The first group was placed at 40 °C (from 25 °C to 40 °C with a 1 °C increment per hour) and maintained for 24 h (2000lx light intensity and 60% humidity). The second group was placed at 40 °C (from 25 °C to 40 °C with a 1 °C increment per hour) and maintained for 24 h (60% humidity). The plants were wrapped with foil to avoid light. The third group was placed at 4 °C (from 25 °C to 40 °C with a 1 °C increment per hour) and maintained for 24 h (2000lx light intensity). The fourth group was placed at 4 °C (from 25 °C to 40 ° C with a 1 °C increment per hour) and maintained for 24 h. (1 degree per hour from 25 °C to 4 °C). The plants were wrapped with foil to avoid light. The drought treatment involved stopping watering for 20 days and measuring the relative soil water content to determine the drought degree (Table S8). Unstressed seedlings were used as the control samples (CK). Three biological replicates were collected from the same position at each collection time for RNA isolation.

RNA extraction and Expression analysis of qRT-PCR

The total RNA extraction was performed using the SPARKeasy A0305 kit (QingDao, China), following the manufacturer’s instructions. Each treatment comprised three replications. The integrity and concentration were detected by 1% agarose gel electrophoresis and a NanoDrop2000 (Thermo Fisher Scientific, Waltham, MA, USA). First-strand cDNA was prepared using the Takara PrimerScriptTMRT reagent kit with gDNA Eraser, and RT-qPCR primers were designed by Primer 5.0 (Table S7). Takara TB Green was used to set up the Premix Ex Taq qGreTli RNaseH Plus for the qRT-PCR experiment. Per the manufacturer’s instructions, 18S was used as an internal reference gene, and there were three biological replicates and three technical replicates in each group. The relative expression was determined by a three-step method (95 °C for 5 s, 55 °C for [15 s], and 72 °C for 30 s, repeated 40 times) on an Applied Biosystems StepOnePlus. The gene expression level was calculated according to the quantitative Ct method. The R package pwr was used to conduct Student’s t-tests. The R packages corrplot (0.84) and circlize (version 0.4.11) were used to conduct correlation analysis.

Identification and characterization of CjCBFs

First, we used the transcriptome data of C. japonica (Naidong) as the local blast resource bank to obtain 11 candidate genes. Further screening incorrectly predicted CBF genes and redundant sequences, and four CBF transcription factors responsive to low temperature were identified. They were named according to their homology with CBF transcription factors in related plants (Camellia sinensis). The nucleotide and amino acid sequences are provided in Table S1. The full lengths of the CjCBF s sequences are 1.82 kb, 1.25 kb, 1.42 kb, and 1.35 kb. The CDSs are between 540 and 837 bp, and the number of amino acids are between 179 and 278. The molecular weights of the four cjCBF proteins are 30282.9 kDa, 19449.25 kDa, 26655.03 kDa, and 26526.23 kDa. The isoelectric point ranges from 4.9 to 9.8. Hydrophilicity ranges from 3.02 to 1.8 (Table 1).

Analysis of conserved motifs

A comparison of the CjCBFs sequences (Fig. S1) revealed that the internal similarity of CjCBFs was 73.96%, and the domain sequence was highly conserved. Comparison of the domain sequences of A. thaliana and C. japonica (Naidong) (Fig. 1A) showed that the domain of C. japonica (Naidong) consists of the AP2 family universal DNA-binding domain (AP2 domain) and two signature sequences (DS-W-L), which were slightly different from the conserved domain of A. thaliana. We then analyzed the conserved motifs of A. thaliana and C. japonica (Naidong) and identified 10 motifs (Fig. 1B) with amino acid lengths between 6 and 50. In particular, Motif10 (WYQGDD) was the shortest (6 residues), while Motif1 (KKVRETRHPIYRGVRQRNSGKWVCEVREPNKKSRIWLGTFPTAEM AARAH) and motif2 (DVAAJALRGRSACLNFADSAWRLPIPESLCPKDIQKAA AEAAEAFRPELC) were the longest (50 residues) and contained the AP2 domain. Only CjCBFs have motif8 (YLDVSGDVAKP), and only A. thaliana has motif7 (TTDHGLDMEETLVE). Only CjCBF1 and CjCBF5 have motif10. The functions of most of the motifs have not yet been elucidated. The similarity of motifs within species indicates that CBFs within species are highly conserved. The composition of different conserved motifs may be related to variation in the functional differentiation of genes.

Analysis of the phylogenetic tree

To characterize the evolutionary relationships of CBF s among species, we constructed a genetic tree (Fig. 2) using 424 CBFs from 19 species, including dicotyledons and monocotyledons. These genes were generally divided into three categories. Monocotyledons and dicotyledons were clearly separated, and dicotyledons were divided into two subclasses: Clade I and Clade II. However, some monocotyledons (e.g., oil palm) were divided into A-II subcluster in Clade I, and rice was also partially distributed in this subcluster. This cluster is close to Clade III, with groups containing only Monocotyledons. This may represent the region of functional overlap of dicotyledon and monocotyledon CBFs. This also indicates that this cluster of CBF ancestral genes had differentiated before the separation of monocotyledons and dicotyledons and thus that the function and pedigree of CBFs became gradually enriched. A-I and B-II were the largest subclusters in Clade I and Clade II, respectively. CjCBFs were divided into A-I subclusters and clustered with A. thaliana. CBFs were highly conserved within species, and CBFs of the same species often clustered together.

GO analysis

GO enrichment analysis of CjCBFs was conducted to better understand the biological function of CjCBFs. CjCBFs were enriched in nucleus (GO:0005634), DNA-binding transcription factor activity (GO:0003700), regulation of transcription, DNA-templated (GO:0006355), transcription regulatory region sequence-specific DNA binding (GO:0000976), glucosinolate metabolic process (GO:0019760), cold acclimation (GO:0009631), and response to cold (GO:0009409). The results of the GO analysis confirmed that CjCBFs play a role in cold acclimation (4 °C) (Fig. 3) (Xie et al., 2010).

Expression model of CjCBFs based on RNA-seq

Based on the data analysis of the full-length transcriptome under treatment with low temperature stress, the standardized FPKM cluster analysis heat map (Fig. 4A) shows that CjCBFs responded strongly to low temperature. Although CjCBF1 played a role in all three gradient treatment groups, it was significantly up-regulated at the 12 h stage, and the strength of the effect was significantly higher in the 12 h stage than in the 24 h and 72 h stages. The response of CjCBF3 to low temperature was basically the same in the three treatment groups. CjCBF2 responded in all three treatment groups but was significantly up-regulated in 24 h and then decreased in 72 h. The expression of CjCBF3 and CjCBF5 was similar to that of CjCBF1, except that CjCBF5 was significantly up-regulated in 24 h but then decreased significantly and remained expressed at a low level. qRT-PCR was used to verify the relative expression of CjCBFs. The pattern was nearly identical, indicating that our transcriptome results were reliable (Fig. 5C). The correlation analysis also supports this finding (Fig. 4B). There is much evidence showing that lncRNA participates in biological processes by regulating gene expression, including responses to stress in plants. Our SMART and next-generation sequencing data were used to predict possible lncRNA transcripts. A total of 147 possible CjCBF-containing transcripts (Fig. 4D) were predicted by lncRNA.

Interactions of CjCBFs

To further study the post-transcriptional regulation mechanism of CjCBFs under cold stress, we characterized the interactions among CjCBF1, CjCBF2, CjCBF3, and CjCBF5 using the Y2Hgold yeast two-hybrid technique. All CjCBFs were completely inhibited when the concentration of Aureobasidin A (AbA) reached 90 ng/µl, while the corresponding AbA concentration for CjCBF2, CjCBF3, and CjCBF5 was 30 ng/µl (Fig. 5A), The results of self-activation verification showed that the self-activation ability of CjCBF1 was significantly higher than that of other CjCBFs, CjCBF2-PGAD, CjCBF3-PGAD, and CjCBF5-PGAD were transformed into yeast Y2HGold cells in CjCBF1-PGBD; PGAD-T and PGBD-lam were transformed into the negative control, and PGBD-T and PGBD-53 were transformed into the positive control. The results of interaction verification showed that the yeast cells co-transformed with CjCBF1 and CjCBF2, CjCBF1 and CjCBF5 could grow normally on the nutrient-deficient medium, indicating that CjCBF1 could interact with CjCBF2 and CjCBF5 (Fig. 5B). To verify the yeast two-hybrid results, we analyzed the subcellular localization of CjCBF1, CjCBF2, and CjCBF5. CjCBF1 and CjCBF5 were located in the nucleus and cell membrane, and CjCBF2 was located in the cell membrane (Fig. 5C). These results were consistent with the yeast two-hybrid results. Several CjCBFs were expressed to aid the ability of Camellia japonica to cope with low temperature stress.

Expression of CjCBFs in five Camellia japonica tissues

To more fully characterize the spatial expression characteristics of CjCBFs in C. japonica (Naidong), the leaves, stems, flowers, buds, and roots of C. japonica (Naidong) were treated with 4 °C for 24 h. The expression of CjCBFs was detected by qRT-PCR. the results showed, In the normal cultured CK group, the expression of 4 CjCBFs was lower, but it was worth noting that the expression of CjCBF1, CjCBF3 in flowers was higher than that of CjCBF2, CjCBF5. This suggests that CjCBF1, CjCBF3 may be involved in flower development. after 24 h of 4 °C treatment, 4 CjCBFs were Significant difference, but only CjCBF1 was highly expressed in all five tissues, especially in the petals, and CjCBF2, CjCBF3, and CjCBF5 were expressed at low levels. The expression of CjCBF2 in roots and flowers was low, and the expression of CjCBF3 and CjCBF5 in stems, leaves, and buds was significantly higher, and the expression of CjCBF5 in buds was higher (Fig. 6). These results suggest that CjCBF1 may play an important role in the response to cold stress in C. japonica (Naidong).

Bioinformatics Analysis of CjCBF1 Cis-Elements

Promoter cis-elements function as transcription factor binding sites and are key for the transcriptional regulation in response to abiotic stress (Moore et al., 2010). As CjCBF1 plays an important role in CjCBFs, we cloned and analyzed the 1,400 bp upstream sequences of CjCBF1. Bioinformatics analysis revealed the sequence to contain 31 cis-acting elements (Fig. 7 and Table S3), including 12 photosensitive responses and related cis-elements, CGTCA and TGACG motifs (involved in methyl jasmonate responsiveness). The sequence was also rich in TC-rich repeats related to stress, and MYB and MYC related cis elements. The results suggest that CjCBF1 has a wide response to abiotic stress and that light may play an important role in the expression of CjCBFs.

Expression of CjCBFs under different stresses and the effect of light

Camellia japonica (Naidong) was treated with high temperature stress, drought stress, low temperature stress, and light stress. Only CjCBF1 and CjCBF5 responded to high temperature (40 °C) stress (Fig. 8A), but the expression of CjCBF1 was significantly lower under 40 °C treatment than under 4 °C for 24 h and decreased slightly under joint high temperature and light stress treatment. Under low temperature and shading stress (Fig. 8B), the expression of CjCBFs decreased significantly compared with the non-shading treatment, which indicated that light may not represent a signal that CjCBFs respond to at high temperature; however, the effect of light on CjCBFs was more important at low temperature. The responses of CjCBF1, CjCBF3, and CjCBF5 to drought stress were weak; only CjCBF2 showed a significant response to drought stress, indicating that the responses of CjCBFs in Camellia japonica to drought stress are not strong (Fig. 8C).