Genome-wide investigation of the AP2/ERF superfamily and their expression under salt stress in Chinese willow (Salix matsudana)

Plant material and salt stress treatment

The salinity stress treatments were carried out on Salix matsudana ‘yanjiang’ and Salix matsudana ‘9901’. Our previous research results showed that ‘yanjiang’ was salt-sensitive variety, while ‘9901’ was salt-tolerant variety (Zhang et al., 2016). The two salix cultivars used in this study were collected from the botany garden of Nantong University (Nantong, China). The two salix cultivars were authorized for only scientific research purpose, and were deposited in school of life science in Nantong University. The stem cuttings (length, 8–10 cm; coarse, 2–3 mm) of two materials were selected for hydroponic rooting for 20 days. The stems with new generated roots were dipped into NaCl solution (150 mM) (only root and part of shoots were immersed in the solution) for 4 h. All root samples were divided into four categories with three biological replicate to do RNA sequencing: Sample1-1/Sample1-2/Sample1-3 (‘yanjiang’ without salt stress treatment, YJ CK), Sample1N-1/Sample1N-2/Sample1N-3 (‘yanjiang’ with salt stress treatment, YJ NT), Sample2-1/Sample2-2/Sample2-3 (‘9901’ without salt stress treatment, 9901 CK), and Sample2N-1/Sample2N-2/Sample2N-3 (‘9901’ with salt stress treatment, 9901 NT).

RNA isolation and Real-time Quantitative PCR (qRT-PCR) analysis

Total RNA was extracted using TaKaRa MiniBEST Plant RNA Extraction Kit (Takara, Dalian, China) from roots according to the manufacturer’s instruction. Four samples (Yanjiang, Yanjiang NT, 9901, 9901 NT) were collected following the same samples treatments procedure as that in RNA sequencing. For each sample, 3 µg of total RNA was used to synthesize first-strand cDNA with SuperScriptII reverse transcriptase (Takara, Dalian, China). For qRT-PCR, the reaction preparation, application parameter settings and quantitative analysis were performed as previously described (Chen et al., 2018). The reactions were performed using the ABI Prism 7000 Real-time PCR system (Applied Biosystems, USA). The Salix purpurea Actin1 gene (SapurV1A.0655s0050.1) were used as reference genes. The gene-specific primers for the 15 selected genes are listed in Table S1.

Genome sequence retrieval

The Populus trichocarpa and Salix purpurea sequences were downloaded from JGI (http://www.phytozome.net/). The Salix matsudana sequences were obtained from our sequencing, and assembly results were obtained by Roche/454 and Illumina/HiSeq-2000 sequencing technologies (Zhang et al., 2020).

Identification of AP2/ERF genes in Salix matsudana and Salix purpurea

The Pfam accession number of AP2 domain is PF00847 (Gathering cut-off value, 20.6). We downloaded the Hidden Markov Model (HMM) profile for the AP2/ERF TFs from the Pfam database (http://pfam.xfam.org/) with Pfam accession number PF00847 as the search keyword. An alternative HMM profile was built by sequence alignment using ClustalW (Larkin et al., 2007). Two HMM profile files were provided as File S1 and File S2. Using an in-house Perl script with two HMM profiles as queries, hmmsearch was carried out by searching the Salix matsudana and Salix purpurea protein databases with default parameters. To validate the putative accuracy of two HMM search results, the candidate protein sequences were checked in three websites: SMART (http://smart.embl.de/#), CDD (https://www.ncbi.nlm.nih.gov/cdd/), and Pfam (http://pfam.xfam.org/). Candidate proteins with positive results from all three websites were selected as AP2/ERF family members of Salix matsudana and Salix purpurea. Additionally, putative AP2/ERF protein characteristics, including length, molecular weight, and isoelectric point, were calculated by the ExPasy site (http://au.expasy.org/tools/pi_tool.html). The genes CDS sequences were listed in File S3.

Phylogenetic analysis and classification of SmAP2/ERF genes

Using an in-house Perl script (domain_xulie.pl), the conserved AP2 core domains of putative SmAP2/ERF proteins were obtained and subjected to multiple sequence alignment using ClustalW (Larkin et al., 2007). To better classify these SmAP2/ERF proteins, 48 AP2 domains from known categories of Arabidopsis and Populus trichocarpa AP2 genes were selected to carry out multiple sequence alignment with SmAP2/ERF proteins, and a phylogenetic tree based on this alignment was built by MEGA 7.0 with the neighbor-joining method with default parameters (Kumar, Stecher & Tamura, 2016). Bootstrap value was set to 1,000. Depending on the phylogenetic tree constructed by SmAP2, PtAP2, and AtAP2 domains, these SmAP2/ERF proteins were classified into different subfamilies and subgroups.

Gene structure and conserved motif structure analysis

The UTR–exon–intron structures of the SmAP2/ERF genes were obtained based on the gene annotation gff3 files we assembled. Using the online website tool Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/), we obtained the gene structure diagrams (Hu et al., 2015).

To characterize the structures of SmAP2/ERF proteins, the online tool MEME (http://meme-suite.org/tools/meme) was used to search for conserved motifs (Bailey et al., 2009). The optimized parameters were employed as follows: any number of repetitions, maximum number of motifs = 10, and the optimum width of each motif was 6–50 residues. The search result file meme.xml was downloaded from the website and opened by TBtools v0.66831 to obtain the gene structure diagram (Chen et al., 2020a; Chen et al., 2020b).

Gene position on chromosomes, and gene tandem and segmental duplication analysis

Using the “Amazing Gene Location from GFF3/GTF File” tool of TBtools, the SmAP2/ERF genes were mapped on 38 chromosomes of Salix matsudana. Because some scaffolds were not assembled onto the chromosomes, not all SmAP2/ERF genes mapped onto the chromosomes (Chen et al., 2020a; Chen et al., 2020b).

Salix matsudana is a tetraploid willow. Gene duplication events are a common phenomenon in the genome. There are two kinds of gene duplications in the genome: tandem duplication events (TDs) and segmental duplication events (SDs). TDs refer to two or more adjacent homologous genes located within 200 Kb on a single chromosome; SDs refer to homologous gene pairs between different chromosomes (Cannon et al., 2004). The gene duplication pairs were identified in TBtools by the “Blast compare 2 Seq [sets] <Big File>” and “Quick McscanX Wrapper” tools. The candidate duplicated genes should have ≥80% coverage and ≥65% similarity. The TDs of SmAP2 genes were revealed on a chromosome using the “Amazing Gene Location from GFF3/GTF File” tool of TBtools. The SDs of SmAP2 genes were visualized by the “Amazing Super Circos” tool of TBtools (Chen et al., 2020a; Chen et al., 2020b).

Divergence time calculation of duplicated genes

After BLASTn analysis of CDS sequences and obtaining duplicated gene pairs, the nonsynonymous substitution rate (Ka) and Synonymous substitution rate (Ks) of gene pairs were calculated by the “Simple Ka/Ks calculator (NG)” tool of TBtools. The divergence time was estimated with the formula: T = Ks/2 λ (Song et al., 2016). The clock-like rate λ value (9. 1 ×10⁻⁹) from Populus was used in the calculation (Lei et al., 2012; Lynch & Conery, 2000).

Collinearity analysis between Salix matsudana and the representative species

To demonstrate the syntenic relationships of the orthologous SmAP2/ERF genes obtained from Salix matsudana and other two selected plants (Populus trichocarpa, and Salix purpurea), the syntenic analysis maps were constructed using the “Amazing Super Circos” tool of TBtools (Chen et al., 2020a; Chen et al., 2020b).

RNA sequencing and a heat map generated by hierarchical clustering

Transcriptome sequencing data of 12 samples were obtained by Illumina HiSeq sequencing. Using TopHat2 software (Kim et al., 2013), the clean reads were mapped to the reference genome sequence of S. matsudana. Gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped (FPKM) (Jin, Wan & Liu, 2017). The FPKM values of all genes from RNA sequencing were available as File S4. Differential expression analysis of two conditions/groups was performed using the DESeq R package (1.10.1). To identify DEGs, fold change ≥ 2 and false discovery rate (FDR) < 0.01 were used as screening criteria. Using the “Amazing HeatMap” tool of TBtools, a graph of the expression level of SmAP2/ERF family genes with hierarchical clustering was generated (Chen et al., 2020a; Chen et al., 2020b).

Identification, phylogenetic analysis, and classification of 364 AP2/ERF TF family members in Salix matsudana

By HMM profile search against the Salix matsudana protein database, a total of 364 full-length AP2/ERF family proteins containing at least one AP2/ERF domain were identified as AP2/ERF superfamily members of Salix matsudana (Fig. 1). The original hmmsearch –domtblout results are listed in Fil S5. The name, protein length, molecular weight, and isoelectric point of individual genes are listed in Table S2.

The phylogenetic relationships of SmAP2/ERF proteins were inferred by multiple sequence alignment of the AP2 domain, which included approximately 50–60 amino acids. The sequence alignment of all AP2/ERF genes showed some conserved amino acids at specific positions, as previously reported (Nakano et al., 2006) (Fig. S1). For example, the WLG element (58th–60th amino acids; 58–60AA) was highly conserved in the ERF and RAV families; alternatively, the conserved sequences from 58–60AA were converted into YLG elements in the AP2 family and HLG element in two sololist members (Liu et al., 2019). In many species, these conserved amino acid profiles contribute to convincing classification of AP2/ERF genes. Based on multiple sequence alignments of 48 AP2/ERF proteins from Arabidopsis and Populus trichocarpa with known categories and 364 Salix matsudana AP2/ERF proteins, we constructed a phylogenetic tree using the neighbor-joining method to explore the phylogenetic relationships of Salix matsudana AP2/ERF proteins. The phylogenetic tree showed that there were 55 AP2/ERF genes that belong to the AP2 family, with 47 genes that encode proteins with two AP2 domains and eight genes (SmAP2-20, SmAP2-25, SmAP2-29, SmAP2-35, SmAP2-36, SmAP2-40, SmAP2-41 and SmAP2-55) that encode proteins with a single AP2 domain (Fig. 1). Additionally, 301 genes that were predicted to encode proteins with a single AP2 domain were members of the ERF family. The ERF family could be further classified into two subfamilies, ERF and DREB. Of the 301 members, 166 and 135 genes belonged to the ERF and DREB subfamilies, respectively. The ERF family genes from Salix matsudana were distributed in B1–B6 subgroups; the DREB family genes from Salix matsudana were classified into A1–A6 subgroups. The gene number and percentage of each subgroup are listed in Fig. 2 and Table S3. Six putative genes were classified as RAV subgroup genes that encode proteins containing one AP2/ERF domain and one B3 domain (Fig. 1). Two genes were designated as Soloist genes, whose AP2/ERF-like domain sequences had lower homology compared with other AP2/ERF genes (Fig. 1).

The AP2/ERF genes number, classification and percentage of different subgroups from five plant species, including the model plant Arabidopsis, Populus, and two Salix plants, are listed in Table S3. The gene name of AP2/ERF genes from Populus trichocarpa and Salix purpurea are listed in Table S4. As a tetraploid plant, the total number (364) of AP2/ERF genes was much larger in Salix matsudana than in the other four species. The number of AP2/ERF genes in Salix matsudana was 2.5-, 1.8-, 1.9-, and 2.1-fold higher than those in A. thaliana (Nakano et al., 2006), Populus trichocarpa (Zhuang et al., 2008), Salix purpurea, and Salix arbutifolia (Rao et al., 2015, respectively. For DREB and ERF subfamilies, the percentage of all AP2/ERF genes in Salix matsudana was similar to those of A. thaliana, Populus trichocarpa, and Salix purpurea, and the percentages of DREB and ERF subfamilies were 38% and 45%, respectively. In Salix arbutifolia, the percentage of DREB (33%) was lower than that of the other four species, whereas the percentage of ERF (50.8%) was higher. In Salix matsudana, the percentage of the AP2 subgroup was highest among all five species (15%) and the numbers of most of gene sub-classifications were doubled, including the Soloist gene; there were two Soloist genes in the Salix matsudana genome. However, no duplications were observed in the RAV subgroup, and only six RAV genes were found in the Salix matsudana genome.

Gene structure and conserved motif analysis

To understand the structural diversity of SmAP2/ERF genes in different clades, a different form of phylogenetic tree of SmAP2/ERF family was constructed and the different subgroups were labled (Fig. 3A). The intron and exon structures of SmAP2/ERF genes were revealed by inputting Gff3 files into TBtools (Fig. 3B). A total of 55 genes of the AP2 subfamily had more exons than ERF and other subfamilies. Apart from three exons in the SmAP2-29 and four exons in the SmAP2-20, other members of the AP2 subfamily contained more than seven exons. The intron number was less than three in many members of the ERF and RAV subfamilies. In total, 215 gene members did not have introns (Fig. 3B). The exon/intron structures of genes that were classified in the same clade were similar. Many gene pairs were found in the phylogenetic tree that potentially originated from allotetraploid evolution of Salix matsudana. Many gene pairs (approximately 70%) maintained the same or similar gene structure during Salix matsudana evolution, which indicated that the SmAP2/ERF genes were conserved at the DNA level after polyploidization.

TF proteins always contain many conserved motifs to activate gene expression. A total of 10 conserved motifs were detected in 364 SmAP2/ERF proteins using the online MEME software, and a block diagram was constructed to characterize SmAP2/ERF protein structure (Fig. 3C, Fig. S2). Motif-4, Motif-1, Motif-2, Motif-3, Motif-5, Motif-7, and Motif-9 were found in the AP2 domain regions. The Motif-5 region covered the region of Motif-4 and Motif-1, whereas Motif-7 included Motif2 and Motif3. Motif-9 is a specific motif that is only found in the second AP2 domain of the AP2 subgroup. Motif-1, Motif-2, Motif-3, and Motif-4 were detected in 90% percent of the ERF subfamily proteins. Thirty proteins of the ERF subfamily lacked one or two motifs of Motif-1–4. An extreme example is SmERF B2-13, which only had Motif-2. Motif-6, Motif-8, and Motif-10 are motifs located outside of the AP2 domain. Motif-6 was primarily found in the AP2 subfamily with only one exception, SmERF B4-4, which was in the ERF-B4 clade. In the AP2 subfamily, members with two AP2 domains had Motif-6 located between the two AP2 domains. Motif-8 was found in 69 proteins of the AP2/ERF family, and its location was adjacent to the carboxyl terminal of Motif-3. Many proteins from the DREB-A1, DREB-A4, DREB-A5 clades had Motif-8. Motif-10 was found in 62 proteins of the AP2/ERF family, with 61 proteins from the ERF subfamily and only one from the AP2 subfamily. Motif-10 was mostly distributed on the proteins from the ERF-B3, DREB-A2, and DREB-A4 clades. The functions of these three motifs need to be elucidated by further experimental analysis.

Besides protein SmAP2-20, the entire AP2 domain was distributed in the amino terminal or in the middle position of the proteins. In the two Soloist genes, only one motif, Motif-2, was found. The conserved motif composition and gene structure of the same subfamily were similar, thus verifying the reliability of the phylogenetic tree classification.

Chromosome distribution and duplication of SmAP2/ERF superfamily genes

The chromosome location of the identified SmAP2/ERF genes was constructed using TBtools. In total, 310 genes from the AP2/ERF superfamily were unevenly distributed on 38 chromosomes (Fig. 4); 54 other genes located on scaffolds were not illustrated in Fig. 4. The chromosome with the largest number of AP2/ERF genes was Chr21, which had 22 genes. Only one AP2/ERF gene each was located on Chr14 and Chr36. On the four chromosomes Chr1, Chr3, Chr22, and Chr27, only two AP2/ERF genes were found. In 38 chromosomes, most of the AP2/ERF genes from different subgroups were arbitrarily distributed, such as five of six RAV genes located on Chr15, Chr37, Chr34, Chr31, and Chr11. Moreover, the two Soloist genes were distributed on Chr29 and Chr5. However, SmERF B3 subgroup members clustered together with three genes as a cluster unit. We found 12 clusters in 12 chromosomes (Fig. 4), which accounted for 62% of the whole SmERF B3 subgroup.

In addition, we also analyzed the tandem duplication events (TDs) of the AP2/ERF genes located within in the 200-kb range of chromosomal regions of the Salix matsudana genome. Eleven TD regions, which included 23 SmAP2/ERF genes, clustered into 11 linkage groups (LGs) of the Salix matsudana genome (Fig. 4). LGs that contained cluster repeat genes were Chr7, Chr8, Chr10, Chr13, Chr17, Chr19, Chr21, Chr24, Chr28, Chr31, and Chr33. All genes of the repeat clusters were localized within a genomic segment of less than 20 Kb; for example, SmDREB A4-20 and SmDREB A4-19 were localized on a 3.6-Kb segment of Chr24. On Chr8, three genes clusters (SmERF B3-6, SmERF B3-7, and SmERF B3-8) located on a less than 12-Kb segment. In 11 tandem repeats, six came from the SmERF B3 subgroup, two came from the SmDREB A4 subgroup, and one each came from the SmDREB A1 and AP2 subgroups. SmERF B3-40 and SmERF B3-39 tandem repeat pairs had 97% protein sequence identity.

In addition to tandem duplications, many segmental duplication events (SDs) were found in Salix matsudana by MCScanX (Fig. 5, Table S5). We found a total of 28,348 collinear gene pairs (not shown) in the Salix matsudana genome, from which 298 AP2/ERF collinear gene pairs were identified. Then, Ka, Ks, and Ka/Ks ratios of these 298 AP2/ERF collinear gene pairs were calculated to estimate the divergence time (T value) and selection pressure among duplicated SmAP2/ERF gene pairs (Table S6). All of the Ka/Ks values were below 1, which indicated that these genes might have experienced strong purifying selective pressure during evolution. Among the 298 AP2/ERF collinear gene pairs, 198 were located on duplicated segments on 38 chromosomes in Salix matsudana (Fig. 5 and Table S3). The collinear gene pairs in the Salix matsudana genome were visualized by Circos, and the gene pairs were linked by lines (grey lines indicated all gene pairs, red lines indicated AP2/ERF collinear gene pairs).

Synteny analysis of AP2/ERF genes between Salix matsudana and two related Salicaceae species, Populus trichocarpa and Salix purpurea

To further infer the phylogenetic mechanisms of the SmAP2/ERF family, we constructed two comparative syntenic maps of Salix matsudana with two related species, Populus trichocarpa (Fig. 6A) and Salix purpurea (Fig. 6B). Collinear AP2/ERF genes pairs between Salix matsudana and other two species are listed in Table S7. A total of 263 SmAP2/ERF genes showed syntenic relationships with 183 genes from Populus trichocarpa, and 248 SmAP2/ERF genes showed syntenic relationships with 144 genes from Salix purpurea. The number of orthologous pairs between Salix matsudana and Populus trichocarpa, and Salix matsudana and Salix purpurea were 423 and 292, respectively (Table S7). Some PtAP2/ERF and SpAP2/ERF genes were found to be associated with at least four syntenic gene pairs. Interestingly, the number of collinear gene pairs identified between Salix matsudana and Salix purpurea were less than that between Salix matsudana and Populus trichocarpa.

In the comparative syntenic map between Salix matsudana and Populus trichocarpa, syntenic links were found between all 19 Populus trichocarpa chromosomes and all 38 Salix matsudana chromosomes (Fig. 6A). Alternatively, in the comparative syntenic map between Salix matsudana and Salix purpurea, there were no syntenic links between Chr1, Chr12, and Chr36 from Salix matsudana, and Chr15Z and Chr15W from Salix purpurea (Fig. 6B).

Specific expression of AP2/ERF superfamily genes under salt stress

To investigate the physiological roles of SmAP2/ERF genes in salt stress tolerance, we identified the expression patterns of SmAP2, SmERF, and SmDREB subgroup genes from the RNA sequencing data. The RNA transcripts of 285 genes were identified from the RNA sequencing data. Using fold change ≥ 2 and false discovery rate (FDR) <0.01 as screening criteria, 68 SmAP2/ERF genes were identified as DEGs. The DEGs names were listed in Table S8. By inputting the FPKM values (Fragments Per Kilobase of transcript per Million fragments mapped) of these genes in TBtools, three heatmaps were constructed using Log₁₀-transformed values of the FPKM values to demonstrate the expression pattern change under salt stress (Fig. 7).

The expression patterns of 285 genes are illustrated in Fig. 7, and included 47 AP2 (Fig. 7A), 108 DREB (Fig. 7B), and 130 ERF subgroup genes (Fig. 7C). In the AP2 subgroup, the Log₁₀-transformed values of 31 genes were <3, which indicated lower expression in the root and no response to salt stress. Five DEGs had differential expression patterns. The expression levels of four genes (SmAP2-38, SmAP2-4, SmAP2-3, and SmAP2-33) were induced by salt stress, whereas the expression of gene SmAP2-15 decreased after salt stress. In the DREB subgroup, 108 genes including 27 DEGs were present in the heatmap, and expression levels of 10 genes, such as SmDREB A1-10, SmDREB A1-9, and SmDREB A1-7, were highly induced by salt stress and remained higher. In the ERF subgroup heatmap, 130 genes (containing 33 DEGs) were included. The expression levels of 13 genes were intensely upregulated by salt stress, but only the expression of SmERF B4-1 was higher. Three genes were downregulated by salt stress, including SmERF B3-52. In many paralog gene pairs, we found one gene with higher expression, whereas the other gene had lower expression, such as SmDREB-A9/SmDREB-A10, SmAP2-33/SmAP2-39 and SmERF-9/SmERF-10 gene pairs. Fourteen genes with upregulated expression patterns were verified by qRT-PCR (Real-time Quantitative PCR) (Fig. 8). From the results, we found that most genes’ expression patterns were consistent with the FRKM values except SmERF B3-42 and SmAP2-15 genes. Both in ‘Yanjiang’ and ‘9901’ samples, the expression level of thirteen genes was induced after salt treatment, but the expression level of seven genes (SmAP2-33, SmDREBA4-24, SmDREBA1-4, SmDREBA1-7, SmDREBA5-23, SmERF B3-45 and SmERF B4-1) was much higher induced in ‘9901’ than that in ‘Yanjiang’, which was not found in RNA sequencing results.