Genome-wide analysis of BpDof genes and the tolerance to drought stress in birch (Betula platyphylla)

Identification of Dof genes from B. platphylla

The assembled birch genome (Chen et al., 2021) was analyzed, and the unigenes were searched using BLASTX against the NR and Swiss-Prot databases for functional annotation (Camacho et al., 2009). All BpDofs (Table S1) were detected from the birch genome by employing a hidden Markov model (HMM) profile of the Dof domain (PF02701) obtained from Pfam (http://pfam.xfam.org) using the HMMER3.0 program (http://hmmer.janelia.org) (Finn et al., 2016). The conserved domains of these putative BpDof proteins were identified by searching against the NCBI’s conserved domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi?) (Marchler-Bauer et al., 2015). Predictions of theoretical molecular size, isoelectric point (pI), and other physicochemical properties were conducted using the ExPASy ProtParam tool (http://www.expasy.org/tools/protparam.html) (Gasteiger et al., 2003). Nuclear localization signals and transmembrane domains were predicted using SubLoc (http://cello.life.nctu.edu.tw/cello2go/) and TMHMM server 2.0 (http://www.cbs.dtu.dk/services/TMHMM/), respectively (Chen, Huang & Sun, 2006; Krogh et al., 2001).

Sequence alignment and phylogenetic analysis of the 26 Dof proteins

The multiple sequence alignments of 26 birch Dof proteins were arrayed in the ClustalW2 of BioEdit7.2.1 software (http://www.ebi.ac.uk/Tools/clustalw2/) (Larkin et al., 2007). A phylogenetic tree of 26 birch Dof proteins with the 39 Dof proteins of Arabidopsis (Table S2) was constructed using MEGA7.0 by the neighbor-joining (NJ) method (Kumar, Stecher & Tamura, 2016; Tamura et al., 2011). The phylogenetic relationships of 26 birch Dof proteins were analyzed.

Gene structure and conserved motif analysis

The genome sequences of the 26 Dof genes were acquired from the birch genome, and their exon/intron structures were analyzed using the GENE Structure Display Serve 2.0 (http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015). Their conserved motifs with default parameters were analyzed using MEME (http://meme-suite.org/tools/meme) (Bailey et al., 2009), but the maximum number of motifs was set as 15.

Chromosomal distribution and synteny analysis

The length of each chromosome and the location of each BpDof gene (Table 1) were retrieved from the birch genome, and the chromosomal distribution of BpDof genes was visualized with the Amazing Super Circos software in TBtools (Chen et al., 2020a). One-step MCScanX SuperFast software in TBtools was used for gene synteny and collinearity analyses with default parameters, and the syntenic map was constructed with the chromosomal locations of BpDofs. Furthermore, to analyze the selection pressure of BpDofs, the non-synonymous rate (Ka), synonymous rate (Ks), and Ka/Ks values of the corresponding BpDofs were calculated using KaKs Calculator v2.0 (Wang et al., 2009).

Cis-element analysis in the promoters of 26 BpDofs

The 2,000-bp length of the upstream DNA sequence with the 5′-untranslated region (UTR) for each Dof gene was obtained from the birch genome. The cis-elements in the promoter sequences of the 26 BpDof genes were analyzed using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

Plant materials and drought stress treatment

Birch seeds were placed in a glass bottle loosely wrapped with gauze and rinsed with water for 3 days. After filtering out most of the water, the seeds were evenly spread on the soil (soil:vermiculite:perlite = 3:1:1). A thin layer of soil was applied to cover the seeds, after which they were watered and sealed using preservative film with some holes in it. The film was removed from the seeds after germination. Once the seedlings grew to 3–4 cm, the seedlings with uniform growth and in good conditions were selected. One birch seedling was planted in one pot. After about 2 months cultivation, the healthy birch seedlings about 25 cm height with similar growth conditions were treated using 20% PEG6000 for 0.5, 1, 3, 6, 12, and 24 h in a reversed time order. The control plants were treated with fresh water for 24 h. Each birch seedling was watered with 1 L of 20% PEG6000 or water. The developing roots, stems, and leaves of the birch seedlings under drought treatments were collected after 24 h. All samples were immediately frozen in liquid nitrogen and stored at −80 °C. Three biological replicates were conducted in each experiment.

RNA isolation and real-time PCR validation

The RNA of each sample was extracted using a Universal Plant RNA Extraction Kit (BioTeke Corporation, China) from 100 mg plant tissues (roots, stems, or leaves) of birch plants, and cDNA was synthesized from approximately 1 µg total RNA using PrimeScript IV First-Strand cDNA Synthesis Mix (TaKaRa, Japan). A 20- µL volume containing 10 µL of SYBR Green Real-time PCR Master Mix (BioTake Corporation, China), 1 µL of cDNA template, and 1 µL each of the forward primer (10 µM) and reverse primer (10 µM) was used (all primers and amplicon sizes are shown in Table S3), after which ultrapure H₂O was used to make up the reaction volume. The amplicons were completed as follows: 94 °C for 30 s; 94 °C for 12 s, 58 °C for 30 s, and 72 °C for 45 s, for 45 cycles, followed by 79 °C for 1 s for plate reading using an qTOWER³ G, analytik Jena AG, Germany. After the final PCR cycle, the temperature of 0.5 °C per second was increased from 55 °C to 99 °C to generate the melting curve for samples. The relative mRNA levels were determined by normalizing the PCR threshold cycle number of each gene to that of ubiquitin (GenBank number: FG065618) and α-tubulin (GenBank number: FG067376) as internal references for all treatments. The threshold for the Ct values was the machine setting, and the average Ct value was calculated using three biological replicates. The relative expression levels of the 26 BpDof genes were calculated from the threshold cycle by the delta–delta CT method (Pfaffl, Horgan & Dempfle, 2002).

Vector construction and transient transformation

The full-length coding sequences (CDSs) of BpDof4, BpDof11, and BpDof17 were amplified by PCR,and then constructed into the pROKII vector digested with Sma I (NEB, USA) using the In-Fusion™ CF Liquid PCR Cloning kit (Takara, Japan) under control of the CaMV 35S promoter for overexpression of BpDof4, BpDof11, and BpDof17 (35S:BpDof), respectively. The primers used for PCR are shown in Supplemental Table S4. The pROKII-35S::BpDof4, 35S::BpDof11, and 35S::BpDof17 were separately transformed into 4-week-old birch seedlings by Agrobacterium tumefaciens–mediated transient expression (Zhang, Wang & Wang, 2012) with some modifications. In brief, Luria-Bertani (LB) liquid medium supplied with 50 mg/L kanamycin and 50 mg/L rifampicin was used to culture the A. tumefaciens strain EHA105 transformed with pROKII-35S::BpDof4, pROKII-35S::BpDof11, pROKII-35S::BpDof17, or the empty pROKII-35S vector. A. tumefaciens cultures were re-suspended in the transformation solution (1/2 MS + sucrose [2.0%, w/v] + 10 mM CaCl₂ + 120 µM acetosyringone + 200 mg/L DTT + Tween-20 [0.02%, v/v], pH 5.8], which were then harvested at an OD₆₀₀ of 0.6 by centrifugating at 3000 g for 10 min. For transient genetic transformation, the plants were soaked into the transformation solution and shaken at 120 rpm and 25 °C for 2 h. Then the plants were planted vertically on 1/2 MS solid medium (1/2 MS + sucrose [2.0%, w/v] + 120 µM acetosyringone + 200 mg/L DTT, pH 5.8) and incubated at 25 °C in the dark. After culturing for 48 h, the plants were assumed to have been transformed and were then used for subsequent experiments.

Stress tolerance analysis of BpDof4-, BpDof11-, and BpDof17- overexpressing plants

The birch plants overexpressing BpDof4, BpDof11, or BpDof17 were treated with 20% PEG6000 for 6 h. The pROKII-35S transformants and the wild-type (WT) birch seedlings were also treated with PEG6000. Water treatment was used as control. The detached leaves of birch plants were incubated with 0.5 mg/mL nitroblue tetrazolium (NBT, dissolved in phosphate buffer, pH 7.8) and 1.0 mg/mL 3 ′-diaminobenzidine (DAB, dissolved in phosphate buffer, pH 3.8) as described previously (Zhang et al., 2011). Evans blue (1.0 mg/mL, dissolved in sterile deionize water) staining was performed to detect cell death following the published protocols (Kim et al., 2003). Superoxide dismutase (SOD) and peroxidase (POD) activities, H₂O₂ content, and electrolyte leakage were measured as previously described (Liu et al., 2015; Wang et al., 2015b). Three independent biological replicates were performed.

Statistical analysis

All statistical analyses were performed using SPSS software (IBM, IL, USA), and ANOVA was used to determine statistically significant differences between results. The level of significance was set at P <0.05.

Identification and characterization of the 26 Dofs in B. platyphylla

Twenty-six full-length Dof TFs (GenBank accession numbers: MW538484 –MW538509) were identified in B. platyphylla (Table 1) using the HMMER3.0 program with a HMM profile of the Dof domain (PF02701), and conserved domains of these putative BpDof proteins were identified by searching against the NCBI’s conserved domain database. These proteins encoded by BpDof genes consisted of 160–562 amino acids (aa). The molecular sizes and pI values of these proteins ranged from 17.8 kDa to 61.0 kDa and 4.97 to 9.36, respectively, and their aliphatic indexes were between 43.34 and 62.44. Analyses of instability indexes and grand average of hydropathicities indicated that most BpDof proteins were unstable hydrophilic proteins except BpDof16 and BpDof20. The charge results showed that BpDof8 and BpDof14 were neutral; BpDof2, BpDof3, BpDof4, BpDof5, BpDof21, BpDof25, and BpDof26 were negative; and the other 17 BpDofs were positive. All BpDof proteins were predicted to be localized to the nucleus. The transmembrane domain analysis indicated that the 26 BpDof proteins did not have α-helical transmembrane motifs (Tables S5 and S6).

Sequence alignment and phylogenetic analysis of BpDof proteins

The multiple sequence alignments of the 26 birch BpDof proteins (Table S1), together with several representative Dof proteins selected from the published Arabidopsis databases (https://www.arabidopsis.org) and Populus trichocarpa v3.0 genomics resource (https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Ptrichocarpa) (Table S7), were analyzed. The single Dof domain with the C(x)2C(x)21C(x)2C zinc finger pattern was harbored in the N-terminal region of all the putative 26 BpDofs (Fig. 1).

The phylogenetic relationships of the 26 BpDof protein sequences, along with the AtDof protein sequences in A. thaliana, were examined using the NJ method. Results showed that the 26 BpDof proteins were classified into four major groups (A to D) (Fig. 2). The subfamily D1 in Group D was the largest group, which included seven members accounting for 26.9% of all BpDof proteins. The subfamily B1 consisted of five members, accounting for 19.2%. The subfamily C2.1 comprised four members, accounting for 15.4%. Group A and subfamilies C1, C2.2, and D2 contained two members accounting for 7.7%, respectively. Subfamilies B2 and C3 with only one member had a proportion of 3.8%, respectively.

Conserved motifs and gene structure analysis

Conserved motifs were identified using the MEME tool, and an unrooted phylogenetic tree was constructed based on BpDof protein sequences (Fig. 3A). Different numbers of conserved motifs were set in MEME so as to find the most significant conserved motifs in the 26 BpDof proteins based on the statistical significance (E value < 10⁻⁵). The results indicated a total of 15 conserved motifs, and the 26 BpDof proteins were classified into four main groups (A–D) including nine subfamilies basing on the phylogenetic tree. Among them, motif 1 and motif 2 were common motifs shared in almost all BpDof proteins except that BpDof20 lacked motif 2, implying that they were conserved motifs. Some of the BpDof proteins possessed specific motifs, for example, motif 7 was only present in BpDof3, BpDof14, BpDof17, BpDof21, and BpDof26, and motif 11 and motif 13 were only distributed in BpDof6, BpDof7, and BpDof22, suggesting these motifs may be relevant to various functions of BpDof genes.

The transcript sequences of 26 BpDofs were compared with genomic sequences to obtain the distribution of introns and exons (Fig. 3B). The number of introns in the 26 BpDofs ranged from 0 to 2. As a result, 12 BpDof genes contained only exons without any introns (BpDof1, BpDof2, BpDof4, BpDof5, BpDof8, BpDof11, BpDof12, BpDof13, BpDof15-BpDof18), 9 BpDof genes (BpDof6, BpDof7, BpDof14, BpDof19-BpDof22, BpDof24 and BpDof25) respectively contained one intron and two exons, whereas 5 BpDof genes (BpDof3, BpDof9, BpDof10, BpDof23 and BpDof26) contained two introns and three exons.

Chromosomal distribution and inter-specific synteny analysis of BpDofs

We mapped the birch Dof family genes on birch chromosomes to obtain their chromosomal distribution. The results (Fig. 4) showed that the 26 BpDofs were unevenly distributed on the 14 chromosomes of birch genome, with only one BpDof gene located on Chr02, Chr03, Chr10, and Chr13, respectively. Two BpDof genes were located on Chr04 and Chr08, respectively. Three BpDof genes were located on Chr07, Chr11, and Chr14, respectively. Four BpDof genes were located on Chr12, and five BpDofs on Chr06.

During plant evolution, gene duplication is a universal event in all organisms and is important in dissecting the novelties (Lynch & Conery, 2000). According to Fig. 4, seven duplication events were predicted among six chromosomes (Chr02, Chr04, Chr06, Chr08, Chr11, and Chr14), and these duplication events occurred in three subfamilies (B1, C2.1, and D1) of BpDof genes. All of the potential duplication events were inter-chromosomal duplication that occurred between two different chromosomes. Furthermore, the duplicated genes belonged to the same subfamilies, and the three groups of genes were found to have strong collinearity. BpChr06G02126 and BpChr02G19918 were in one group, BpChr06G09621, BpChr11G09292, and BpChr14G12625 were in another group, and the last group included BpChr04G00494, BpChr08G17028, and BpChr11G05806. The Ka/Ks values of these BpDof genes were all less than 1 except those of BpChr11G09292-BpChr14G12625 (Table 2), indicating that they have undergone strong purifying selection during the evolution. The Ks value of BpChr11G09292-BpChr14G12625 was NaN leading to Ka/Ks value as NaN, indicating that gene duplication caused mutation at the nucleic acid level but not at the amino acid level.

Analysis of promoter cis-elements of BpDofs

The putative cis-elements within the 2,000-bp genomic sequences upstream with the 5′- UTR of 26 BpDofs were predicted in the PlantCARE database. The results indicated that some cis-elements in 26 BpDof gene promoter regions were identified, including MYB and MBS (MYB-binding site, involved in drought inducibility), TC-rich repeats (defense and stress-responsive element), LTR (low temperature-responsive element), GT1-motif (light-responsive element), G-Box (light responsiveness element), GA-motif (light-responsive element), W-Box (WRKY-binding site, involved in abiotic stress responsiveness), ABRE (abscisic acid–responsive element), ERE, and GARE-motif (gibberellin-responsive element) (Fig. 5).

Expression patterns of BpDof genes in response to drought stresses

The 2-month-old uniform seedlings were subjected to stress treatment to explore the change in birch Dof expression levels under drought stress. The relative expression levels of the 26 BpDof genes were analyzed in the roots, stems, and leaves of birch under 20% PEG6000 treatment compared with the control (water treatment) (Fig. 6 and Table S8).

In roots, the 26 BpDof genes were differentially expressed under PEG6000-simulated drought stress; most of them were upregulated at nearly all time points, and only a few genes at several time points were slightly downregulated. Compared with the control, the expression of BpDof11 and BpDof17 was significantly upregulated by about 4-fold to 199-fold, and about half of BpDofs reached their peak expression at 12 h after treatment. In stems, the expression of most BpDofs showed no obvious change. The expression of BpDof3 (6-fold) at 12 h, BpDof22 (7-fold) at 12 h, and BpDof24 (7-fold) at 6 h was slightly upregulated. However, the expression of BpDof5 was largely downregulated by about 26-fold at 12 h compared with the control. In leaves, the expression of most genes was significantly different compared with the control at most treatment time points. For example, the expression of BpDof4 (322-fold), BpDof5 (68-fold), and BpDof14 (106-fold) peaked at 0.5 h, while that of BpDof17 (59-fold) peaked at 24 h. However, the expression of BpDof24 was downregulated at all treatment time points compared with the control.

Plants overexpressing BpDof4, BpDof11, and BpDof17 had reduced oxidative stress and cell membrane damage

The expression of BpDof4, BpDof11, and BpDof17 was markedly upregulated at some times under drought treatment compared with the control. Therefore, these three genes were selected for further exploration. NBT and DAB staining was used to detect O²⁻ and H₂O₂ levels so as to determine whether the drought tolerance was strengthened in the BpDof4-, BpDof11-, and BpDof17 - overexpressing plants (Fig. 7). The leaves from BpDof4, BpDof11, BpDof17, WT and pROKII-35S plants were stained with NBT or DAB; the stained leaves from water-treated plants were used as controls. Under drought treatment, the O²⁻ and H₂O₂ levels in the leaves of BpDof4-, BpDof11-, and BpDof17- overexpressing plants were greatly reduced compared with those in WT and pROKII-35S plants. The contents of O ²⁻ and H₂O₂ negatively reflected the reactive oxygen species (ROS) scavenging ability of plants. Therefore, the results indicated that BpDof4-, BpDof11-, and BpDof17-overexpressing plants had enhanced abilities to scavenge ROS. In addition, Evans blue staining was used to detect cell membrane damage. Further, BpDof4-, BpDof11-, and BpDof17-overexpressing plants showed less intense blue staining compared with WT and pROKII-35S plants under drought treatment, suggesting that they had decreased cell death.

Physiological characterization of BpDof4, BpDof11, and BpDof17-overexpressing plants

In this study, SOD and POD activities, H₂O₂ content, and electrolyte leakage were used to assess the resistance of BpDof4-, BpDof11-, and BpDof17 - overexpressing plants to drought stress, as well as that of the pROKII-35S transformants and WT birch plants (Fig. 8 and Table S9). Water treatment was used as control. SOD and POD play important roles in removing ROS in plants under stress. Results showed that the SOD activity in BpDof11 - overexpressing plants was significantly higher than that in WT and pROKII-35S plants; however, the SOD activity in BpDof4- and BpDof17- overexpressing plants had not obvious change compared with WT and pROKII-35S plants. POD activity was significantly higher in BpDof4-, BpDof11-, and BpDof17- overexpressing plants compared with WT and pROKII-35S plants under drought stress. The H₂O₂ level decreased by 19.13%, 36.12%, and 7.01%, respectively, in BpDof4-, BpDof11-, and BpDof17 - overexpressing plants compared with the control. Cell death was evaluated using the electrolyte leakage rate. Electrolyte leakage showed a slight decrease in BpDof4-, BpDof11-, and BpDof17 - overexpressing plants compared with the pROKII-35S and WT plants under drought stress. These results suggested that BpDof4, BpDof11, and BpDof17 could enhance the ROS scavenging ability and inhibit cell death in plants.