Structure and expression analysis of seven salt-related ERF genes of Populus

Screening and phylogenetic analysis

In recent years, we have focused exclusively on the expression and function of transcription factors of Populus in response to abiotic stresses (Wang et al., 2014a; Yao et al., 2016b; Wang et al., 2018). Seven ERFs induced by multiple abiotic stresses were selected for the detection of expression patterns related to salt stress in this study (Yao et al., 2019). The systematic search of the seven ERFs in Populus was performed using PlantTFDB (Version 4.0, http://planttfdb.gao-lab.org/) (Jin et al., 2017). The orthologous sequences in Arabidopsis were searched and downloaded from the Arabidopsis genome TAIR 9.0 (http://www.Arabidopsis.org/index.jsp). Multiple sequence alignments of ERF proteins were performed using Clustal X 1.83 and BioEdit 7.0.5.3 (Thompson et al., 1997; Hall, 1999). Phylogenetic trees constructed with MEGA 10.0.5 using the Neighbor Joining (NJ) method were carried out with 1,000 iterations bootstrap test (Kumar et al., 2018). The information of these seven ERF TFs is listed in Table S1.

Gene structure and conserved motifs

Gene structure of seven ERF TFs were detected using TBtools (http://www.tbtools.com/) by comparison with their corresponding genomic DNA sequences from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Ptrichocarpa) (Chen et al., 2018). Conserved motifs in seven Populus ERF TFs were detected using the MEME software, version 5.0.2 (Bailey et al., 2009), which was run using parameters: any number of repetitions; maximum number of motifs, 6; optimum motif width, between 6 and 50 residues.

Chromosome location and synteny analysis

The chromosome location of these seven TFs were illustrated using the “Map Genes On Genome From Sequence Files” methods of TBtools (Chen et al., 2018). The subject sequences used in this study were downloaded from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Ptrichocarpa). Gene synteny of intra-species and collinearity characteristics of the seven TFs were analyzed using MCScanX with the default parameters and considering pBLAST ≤ 1e⁻⁵ (Wang et al., 2012). Colinear blocks between two sets of linkage groups (LGs) are linked using Circos-plots tool of TBtools (Chen et al., 2018). In this study, all of the 209 ERFs of Populus were obtained from PlantTFDB and used for synteny analysis (Jin et al., 2017).

Cis-elements of promoter

The 2,000 bp upstream sequence of the translation start site were downloaded from the Phytozome v12.1 database and signed as the promoter sequence, respectively. The prediction and determination of the location of cis-elements in promoter sequence was performed with the PlantCRAE software (Lescot et al., 2002).

GO annotation and enrichment

GO annotation and enrichment of seven ERF TFs were performed with OmicsBox 1.2.4 (https://www.biobam.com/omicsbox). According to the manual, blast annotation was proceed with the default configurations. On the basis of GO classification, these seven ERF TFs were involved into biological processes, molecular functions and cellular components three terms. GO enrichment was carried out using Fisher’s exact test methods in this study.

Plant materials and salt stress treatment

The hybrid aspen 84K poplar (P. alba × P. glandulosa) were used as the materials in this study. Twigs from the same clone were cut and planted in Murashige and Skoog (MS) medium for 20 days under the controlled conditions: relative humidity, 60%–70%; light/dark cycle, 16/8 h; average temperature, 25 ± 2 °C. To monitor changes in gene expression patterns according to the different treatments, a time-course experiment was designed in this study. Strong and healthy strains with a similar state were subjected to the following treatments: 0.05 M NaCl for 0, 2, 12, 24, and 48 h. The first time point (hour 0) served as a control. After the completion of each treatment, root, stem, leaf, and shoot tissues were harvested from six seedlings, respectively. The samples from six seedlings per time point were pooled per tissue type, frozen immediately in liquid nitrogen, and stored at −70 °C for RNA isolation and RT–qPCR analysis. Total RNA was extracted according to the manuals (Wang et al., 2018).

RT–qPCR

RT–qPCR was performed as described in our previous study (Wang et al., 2014a). Two housekeeping genes, Actin and EF1, were used as reference to monitor the expression change of target genes (Regier & Frey, 2010). The primers used for RT–qPCR are listed in Table S2. The relative expression level was calculated as 2^−ΔΔCt.

Validation of the expression data using RNA-Seq

RNA sequencing was used to cross-validate the gene expression data of RT–qPCR and the results of the two methods were compared. Twenty-six-day-old seedlings with the same genetic background and similar state of the 84K poplar and grown in MS medium were divided into two groups and treated with 100 mM NaCl (S1, S2, and S3) or regular water (W1, W2, and W3) for 24 h, respectively. Roots, stems, leaves, and shoots were respectively collected from each of the replicates. Three biological replicates were tested. Total 24 samples were collected and shipped to Majorbio Bio-pharm Technology Co.,Ltd. (Shanghai, China) for RNA-Seq using the Illumina Novaseq 6000 platform. Sequencing-library construction and RNA-Seq data analysis were proceed referred to previous study (Wang et al., 2014a). The results of gene expression level are indicated as transcripts per million reads (TPM).

Statistics

A single-variable analysis was used to compare gene expression under salt and normal conditions, as well as between different tissues, using a t-test. We used an unsupervised clustering analysis to identify genes with similar expression patterns based on RNA-seq data. A correlation analysis was used to assess the relationships of genes. The data represent the mean ± SD of expression level. Statistics was calculated in R (v3.5.1, http://cran.r-project.org/).

Screening and phylogenetic analysis

In this study, seven ERF genes of Populus that were reported to be related to abiotic stresses (Yao et al., 2019) were selected as the targets for the detection of expression patterns under salt stress. These genes were designated as PtERF001–PtERF007, for convenience (Table S1). The average length of the proteins encoded by them was about 215 amino acids, with molecular weights ranging from 14,649.4 to 54,060.9 Da and isoelectric point values ranging from 4.7314 to 10.2631 (Tables S1). Orthologs of the seven PtERFs were searched in Arabidopsis, yielding five hits with identification percent about 35%–64% (Tables S1). PtERF002 and PtERF003 hit the same ortholog, AtERF1, in Arabidopsis, whereas PtERF004 and PtERF005 hit AtERF016 (Table S3). Detailed information on the seven PtERFs is provided in Tables S1 and S3, including TF_ID and orthologs in Arabidopsis, cDNA and protein sequences.

Phylogenetic relationships between these seven PtERFs and the five Arabidopsis orthologs were detected via multiple alignments of the full-length protein sequence using NJ algorithm methods. The seven PtERFs were divided into three subfamilies: ERF-a, ERF-b, and ERF-c (Fig. 1A), with three PtERFs each belonging to subfamilies ERF-a and ERF-c, respectively. Only one TF, PtERF007, was attributed to subfamily ERF-b. Multiple alignment conducted using the AP2 domain sequence showed similar profile with that of the full-length amino acid sequences (Fig. 1C). On the other hand, gene structure and expression patterns of these seven PtERFs supported the classification of the subfamily in this study.

Gene structure and conserved motifs

The structures, including the UTRs, CDS (coding sequence) and introns, of seven ERFs were investigated. Genes in the same subfamily shared similar exon/intron structures in terms of number of introns, with the exception of PtERF003 from subfamily ERF-a, which possessed one intron (Fig. 1B). In addition, PtERF007 in subfamily ERF-b also harbored one intron. The conserved motifs that were predicted using MEME showed that motif-1 was conserved among all seven PtERF proteins (Fig. 1B). Motif-3 was only found in members of subfamily ERF-a, whereas motif-5 was only contained in subfamily ERF-c. PtERF001 of subfamily ERF-c did not have motif-5, but carried motif-6. These findings indicate that the structure and function of the poplar ERF TFs are similar, especially among those belonging to the same subfamily.

Analysis of chromosome location and synteny of ERF genes

In this study, seven PtERF genes were mapped on seven out of 19 LGs (Fig. 2). A synteny analysis of a total of 73013 genes of Populus was performed using MCScanX, which led to the identification of 21519 (29.47%) collinear genes (Fig. 2, Table S4). In contrast, 83 (39.71%) collinear genes were identified among the 209 ERF genes and considered to be the result of a whole-genome duplication event. They were located mainly on eight LGs of Populus (Fig. 2, Table S6). Our target gene PtERF001 was included in these 83 collinear genes. In this study, 30 ERF transcription factor genes of poplar, including PtERF002, PtERF004, and PtERF006, related to tandem duplications that impact the expansion of the ERF gene family were identified (Fig. 2, Table S6) and were distributed on 12 of the 19 LGs (Table S5).

Cis-elements in promoter sequences

Putative cis-elements in the promoter region of the seven PtERFs were detected using PlantCRAE. The core-sequence and function annotations of each putative cis-element are listed in Table S7. Most of the promoters of these target ERFs had multiple stress-related cis-acting regulatory elements, such as the ABRE, LTR, MBS, and TGACG motifs; however, the type and number of cis-elements varied among the different genes (Table S7). For example, the ABRE cis-element existed in all but PtERF007, and the CGTCA motif was not found in the promoter region of PtERF001. Whereas, the number of the same cis-element of promoter varied among these seven Populus ERF genes (Dataset S1).

GO annotation and enrichment

GO annotation results showed that these seven ERF genes functioned mainly in the cellular_component and molecular_function terms (Fig. 3A and Fig. S1). None of the genes were predicted to play a role in the biological_process term. The GO enrichment analysis showed that these seven ERF genes were highly enriched in heterocyclic compound binding and organic cyclic compound binding activities, although they were also enriched in the DNA-binding transcription factor activity and transcription regulator activity, which may be involved in the plant stress response process (Fig. 3B, Table S8). Moreover, cellular_component prediction showed these seven ERF genes function mainly in intracellular organelles (Fig. 3, Table S8).

Expression patterns under salt stress

The spatio-temporal expression patterns of the seven ERFs were detected using RT–qPCR. The expression of PtERF003, PtERF004, and PtERF006 in roots treated with salt stress was suppressed to varying degrees (Figs. 4C, 4D, 4F), while PtERF002 was overexpressed (Fig. 4B). The expression of PtERF001 and PtERF005 showed a dynamic pattern of “suppressed–induced amplification”, although the suppressed/induced degree and initial time points were different (Figs. 4A, 4E). However, PtERF007 displayed the opposite expression pattern compared with PtERF001 and PtERF005, which was induced at the 2 h time point and was suppressed thereafter (Fig. 4G). In stems, all seven ERF genes were induced by salt stress across the experimental period, with the exception of PtERF004, which was only induced at the initial and final time points (Fig. 4). In leaves, the expression of PtERF002 and PtERF005 was increased under salt stress, while PtERF003, PtERF004, and PtERF006 were downregulated. The expression patterns of PtERF001 and PtERF007 exhibited similarities in that they were overexpressed at the 48 h time point, but PtERF007 was suppressed at the initial time point (2 h). In shoots, the variation in the expression of these seven ERF genes was more complex compared with the tissues described above. The expression levels of PtERF003 and PtERF005 increased after treatment with salt stress, whereas the expression of PtERF004 decreased significantly (Table S9). The expression of PtERF001 and PtERF007 showed an “induced–suppressed” profile and PtERF006 displayed an “induced–suppressed–recovered” pattern; in contrast, the expression of PtERF002 was relatively unstable. In general, these seven ERF genes showed specific, albeit similar, expression patterns.

The tissue-specific expression pattern of the seven ERFs was also investigated in this study. The expression level of each gene in the root tissues under normal conditions (0 h) and salt stress (24 h) conditions was used as the control and normalized as 1.0, respectively (Fig. 5). Under normal conditions, the tissue-specific expression patterns of PtERF001 and PtERF002 were similar, with lowest expression in stems and highest expression in shoots. The remaining five genes were significantly overexpressed in roots (Fig. 5, Table S10). After treatment with salt stress for about 24 h, the relative expression levels of individual genes in different tissues varied significantly, with the exception of PtERF001 and PtERF005 (Fig. 5, Table S10).

Validation of the expression using RNA-Seq

The expression of these seven ERF genes of Populus was also cross-validated using RNA-seq in this study. In root tissues, the expression (TPM) of PtERF001, PtERF002, PtERF004, and PtERF005 was induced after exposure to salt stress, while PtERF003 and PtERF007 were suppressed and PtERF006 showed no difference compared with the control (Fig. 6, Table S11). In stems, the expression patterns of these seven genes were respectively similar to those observed in the roots, with the exception of PtERF003, PtERF006 and PtERF007. It is worth noting that all seven ERF genes were overexpressed in stems and leaves after treatment with salt stress; the expression of PtERF007 in stems and PtERF006 in leaves in particular was about 3.8 times and 4.4 times that of the control under normal conditions, respectively (Fig. 6, Table S11). In shoot tissues, PtERF006 and PtERF007 were significantly induced by salt stress, whereas PtERF004 was repressed by this treatment (Fig. 6, Table S11). The remaining genes showed no significant differences in expression levels between the tested groups in shoots. Tissue-specific expression patterns were also identified using RNA-seq data. The tissues with high gene expression remained high after the stress treatment (Fig. 7A). In addition, the seven ERF genes could be divided into two subgroups based on their expression patterns: PtERF003, PtERF004, and PtERF005 belonged to one subgroup, while the remaining four genes belonged to the other subgroup (Fig. 7A). The expression correlations detected among these seven genes showed that PtERF004 was correlated with PtERF005, while PtERF002, PtERF003, and PtERF007 were co-expressed significantly (Fig. 7B, Table S12).