Genome-wide identification and expression analysis of the ERF transcription factor family in pineapple (Ananas comosus (L.) Merr.)

Genome-wide analysis of the ERF family in pineapple

The AP2/ERF amino acid sequences of pineapple and Arabidopsis were downloaded from Phytozome12 (https://phytozome.jgi.doe.gov/pz/portal.html) (Goodstein et al., 2012). The Hidden Markov Model (HMM) profiles of the AP2 domain (PF00847) and the B3 domain (PF02362) were obtained from PFAM database (http://pfam.sanger.ac.uk/) (Finn et al., 2008) as the queries for search predicted AP2/ERF proteins in the pineapple dataset using HMMER software 3.0 (http://hmmer.wustl.edu/) with a threshold of e-value<e⁻⁵ (Eddy, 2011). BLAST searches were used to determine the predicted AP2/ERFs in pineapple database with all the Arabidopsis AP2/ERFs as queries. All of the candidate genes were further examined by Simple Modular Architecture Research Tool (SMART, http://smart.embl-heidelberg.de/) (Letunic & Bork, 2018) to confirm the protein sequences derived from the selected pineapple AP2/ERF domains. ExPasy proteomics server (http://expasy.org/) (Gasteiger et al., 2003) was used to predicate the basic information of AcoERF proteins, including molecular weight (MW) and isoelectric points (pI).

Gene ontology (GO) annotation

To investigate the putative function of the candidate AcoERF genes, the gene ontology (GO) annotations of the pineapple ERF family members were downloaded from Phytozome12 (https://phytozome.jgi.doe.gov/pz/portal.html) (Goodstein et al., 2012).

Sequence alignment and phylogenetic analysis

Multiple sequence alignment of the AcoERF proteins was performed using DNAMAN (version 9) with the default parameters (Wang, 2015), and the diagram was visualized using ESPript 3.0 (http://espript.ibcp.fr/ESPript/cgi-bin/ESPript.cgi) (Gouet, Robert & Courcelle, 2003). To explore the evolution relationships of ERF gene family members in pineapple, the ERF amino acid sequences from pineapple and Arabidopsis were used. The multiple sequence alignments of all ERF proteins were performed by using MUSCLE (http://www.ebi.ac.uk/Tools/msa/muscle/) (Edgar, 2004), with the default parameters. Subsequently, the phylogenetic trees were constructed using MEGA6.0 software (http://www.megasoftware.net) (Tamura et al., 2013) via the Neighbor-joining (NJ) method with the following parameters: node robustness was detected using the bootstrap method, and the bootstrap was set to 1000 replications. Finally, the phylogenetic trees were generated by iTOL (https://itol.embl.de) (Letunic & Bork, 2016).

Gene structure and conserved motif predictions

The exon/intron structures of AcoERF genes were obtained from the online program Gene Structure Display Server (GSDS, http://gsds.gao-lab.org/) (Hu et al., 2015). The conserved motifs in the AcoERF proteins sequences were investigated by Multiple EM for Motif Elicitation (MEME, http://meme-suite.org/tools/meme) (Bailey et al., 2015). The optimized parameters were as follows: maximum number of motifs was set to 20, the optimum width of motifs was set to 10-50 residues, and other options were set to default.

Chromosomal localization and synteny analysis

The information of chromosome localization of all pineapple ERF genes was obtained from Phytozome12 (https://phytozome.jgi.doe.gov/pz/portal.html). And then, each genes was mapped on chromosomes using MapChart software based on positional information in the pineapple genome project (Voorrips, 2002). For synteny analysis, the sequence of rice and Arabidopsis ERF proteins were downloaded from Phytozome12. The potential anchors within pineapple, between pineapple and Arabidopsis or rice were searched by BLASTP with the score value of E<1e⁻⁵ and top 5 matches (Zhang et al., 2018). The syntenic block of ERFs within pineapple was visualized using the Circos Program (version 0.69) (http://circos.ca/) (Krzywinski et al., 2009).

The syntenic blocks between pineapple and 2 representative plants were constructed using J. Craig Venter Institute (JCVI, https://github.com/tanghaibao/jcvi) (Tanenbaum et al., 2010). The Synonymous (Ks) and non-synonymous (Ka) values of AcoERF sequences were calculated using bio-pipeline-master (https://github.com/tanghaibao/bio-pipeline) (Tang et al., 2008).

Expression analysis of AcoERF genes in different tissues

To investigate the expression patterns of the AcoERF genes in different development and growth processes of pineapple, the RNA-Seq data of pineapple roots, leaves, flowers, fruits (stage 1-6) were downloaded from the online website Pineapple Genomics Database (PGD, http://pineapple.angiosperms.org/pineapple/html/download.html) (Ming et al., 2015). The criterion of different stage of ovule, sepal, stamen and petal was referenced to the previously described methods (Su et al., 2017). Total RNA of these tissues was extracted using RNA extracted kit (Omega Bio-Tek, Shanghai, China) based on the supplier’s instruction, and the samples were prepared based on the published method (Chen et al., 2017). The cDNA libraries for sequencing were constructed using the NEBNext Ultra™ RNA Library Prep Kit for Illumina (NEB). The reads were aligned to pineapple genome by the TopHat 2.0.0 software with default parameters (Trapnell et al., 2012). The Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) value of each gene was calculated using Cuffdiff and Cufflinks software with default parameters. Finally, a heatmap of AcoERF expression profile was produced by the heatmaply R package (Galili et al., 2018).

Analysis of cis-elements from AcoERF promoters

The promoter sequences of the AcoERF genes were downloaded from the Pineapple Genomics Database (http://pineapple.angiosperms.org/pineapple/html/index.html). And then, the stress- related cis-elements of AcoERF genes were predicted using the Plant Cis-Acting Regulatory Element (PlantCARE, http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

Plant material and abiotic treatments

The pineapple materials were provided by the Qin Lab (Center for Genomics and Biotechnology, Fujian Agriculture and Forestry University, Fujian, China). The seedings were grown in a greenhouse at 25 °C, 16-h light/8-h dark photoperiod and a relative humidity of 70% for one month. And then the one-month-old pineapple plants were exposed to the following treatments: cold stress (4 °C), drought stress (350 mM Mannitol) and salt stress (150 mM NaCl). The samples of pineapple roots and leaves were collected after 6, 12, 24 and 48 h treatments. The one-month-old plants that without any stress treatment were used as controls in the expriment. The leaves and roots were harvested at the indicated time, and stored at −80 °C for subsequent analysis.

RNA isolation and qRT-PCR analysis

Total RNA of the pineapple sample was extracted from roots and leaves using plant RNA extraction kit (OMEGA, Shanghai, China) according to the manufacturer’s protocol. And then, total RNA was reverse-transcribed using AMV reverse transcriptase (Takara, Japan) based on the supplier’s instruction. qRT-PCR analysis was conducted using SYBR Premix Ex Taq II system (Takara, Japan) and Bio-Rad Real-time PCR system (Foster city, CA, USA). The reactions were performed in a 20 µl volume containing 10 µl of 2 × SYBR Premix, 8.2 µl of RNase free water, 1 µl of template, 0.4 µl of each specific primer. The primers used for RT-PCR were showed in Table S1 . The RT-PCR reactions program was completed with the following conditions: 95 °C for 30 s; 40 cycles of 95 °C for 5 s and 60 °C for 34 s; 95 °C for 15 s (Cai et al., 2017; Cai et al., 2019). The analyses were confirmed in triplicate. The relative expression level of each AcoERF gene was calculated based on the comparison threshold period (2^−ΔΔCt) method (Pan et al., 2012). The protein phosphatase 2A gene from pineapple was used as the reference gene (Sang et al., 2013). The final data were subjected to an analysis of variance test.

Genome-wide identification and characterization of pineapple ERF family

A total of 103 pineapple genes were identified as possible encoding proteins containing AP2/ERF domain (Fig. 1A). Among these, 27 genes containing two conserved AP2/ERF domains and 2 genes possessed a single AP2/ERF domain together with a B3 domain were assigned to AP2 and RAV family, respectively; Remaining 74 genes with a single AP2/ERF domain were grouped into ERF family, including 20 DREB and 54 ERF subfamily members. Earlier, the DREB subfamily members of pineapple had been divided into 5 groups: I, II, III, IV and V (Chai et al., 2020). To maintain uniformity, all the pineapple ERF genes were provisionally named as AcoERF1- AcoERF74 according to the gene ID in ascending order (Table S2). The full-length of the 74 AcoERF proteins range from 105 (AcoERF36) to 605 (AcoERF17) amino acid residues, while the coding sequence (CDS) length ranges from 318 to 1819 bp, and relative molecular weights (MW) range from 11.21 to 66.74 kDa. The predicted isoelectric points (pI) varied from 4.71 (AcoERF18 and AcoERF52) to 9.63 (AcoERF30) (Table S2). Moreover, the Go annotation of AcoERFs included biological process (GO: 0006355) and molecular function (GO: 0003700) (Table S2). Within biological process, 72 AcoERF genes were associated with regulation of transcription and DNA-templated. For the molecular function categories, 73 AcoERF genes were involved in transcription factor activity and sequence-specific DNA binding.

Multiple sequence alignment and phylogenetic analysis of pineapple ERF genes

Multiple sequence alignment of the 74 AcoERFs indicated that most ERF family members possessed conserved YRG and RAYD elements within the AP2/ERF domain region (Fig. S1). ERF subfamily proteins containing Ala at position 14 and Asp at position19. However, the residues at position 14 and 19 of DREB subfamily proteins is Val and Glu, respectively. Evolution of ERF TFs was further explored based on the phylogenetic tree constructed with pineapple and Arabidopsis ERF genes (Fig. 1B). Earlier, seven Arabidopsis ERF proteins having low homology with other member of ERF family, were assigned to group VI-L and Xb-L according to the analysis of gene structure and conserved motif (Nakano et al., 2006). In this study, these seven Arabidopsis ERF proteins were also used for phylogenetic reconstruction and formed group G and L. The Arabidopsis ERF proteins of group I and III were from the same group in previous study (Nakano et al., 2006). Therefore, group I and III in this study were renamed as C1 and C2, respectively. The result showed that the AcoERFs were divided into 13 groups, except for AcoERF21. Group I contained the most AcoERFs, accounting for 16.22% of total ERF genes. Followed by group J, which has 10 AcoERFs. However, group L has only one AcoERF.

Gene structure and motif composition of pineapple ERF family

To gain more insight into the evolution and structural diversity of the ERF family in pineapple, we analyzed the gene structure and conserved motif of the AcoERF genes. As shown in Fig. 2, 49 AcoERF genes have no intron, accounting for 66.22% of the total number of AcoERF genes, and all genes of group J have no intron. Moreover, the number of introns in these genes varied from 1 to 16, and most AcoERF genes have 1 to 3 introns, whereas AcoERF17 has the maximum number of introns. All genes of group H have 2 introns, except for AcoERF8. To further explore the divergence and functional relationship of AcoERF proteins, a total of 10 conserved motifs in the pineapple ERFs were identified by the MEME software, and the height of each letter in the logo was proportional to the conservation level of amino acid in all sequences analyzed (Fig. S2). As displayed in Fig. 3, motif 1 is present in every pineapple ERF protein, and almost all the proteins contain motif 2, except for AcoERF38 and AcoERF69. In addition, the motifs in different groups indicated that varying degrees of divergence among them. For example, motif 8 is unique to group H. Motif 6 is only present in group H and I. Motif 10 is only present in the members of group C1 and C2. In general, the pineapple ERF proteins in the same group usually contained similar motifs, which indicates that they may play similar roles in the development and growth of pineapple.

Chromosomal distribution and synteny analysis of AcoERF genes

According to Fig. S3, 71 AcoERF genes were unevenly distributed on the 23 pineapple linkage groups (LG). Among them, LG02 contains the most diverse AcoERF genes, accounting for 12.16% of total genes, whereas LG04, LG14, LG19, and LG25 possessed only one AcoERF gene, respectively. Additionally, AcoERF71, AcoERF73 and AcoERF74 exist on unanchored scaffolds, so they could not be conclusively charted on any pineapple linkage groups (Table S2).

Segmental duplication provide an important mechanism for the expansion of gene families (Lynch & Conery, 2000; Vision, Brown & Tanksley, 2000; Cannon et al., 2004; Zhang et al., 2018). A total of 16 segmental duplication evens with 24 AcoERF genes were found from pineapple genome, such as AcoERF13/AcoERF5, AcoERF16/AcoERF31, AcoERF19/AcoERF35, all of these were segmental duplicates (Fig. 4, Table S3), and these segmental duplication pairs play crucial roles in the expansion of ERF family in pineapple genome. In order to further study the phylogenetic and evolutionary trait of ERF family in pineapple, we set up 2 comparative syntenic maps of pineapple associated with 2 representative plant, including Arabidopsis and rice (Figs. 5A, 5B). According to the syntenic results, there are 37 and 44 orthologous gene pairs in Arabidopsis and rice, respectively (Tables S4 and S5). Some AcoERF genes had multiple orthologous gene pairs (one pineapple gene associated with multiple Arabidopsis or rice genes), such as AcoERF62 associated with AtERF008/019/020, and AcoERF29 associated with OsERF025/026/027, suggesting these genes might play crucial roles in the evolution of ERF family. Some AcoERF genes could be found orthologous gene pairs in Arabidopsis and rice, such as AcoERF15, AcoERF22 and AcoERF32, suggesting that these genes might already exist before speciation. Some AcoERF genes only existed orthologous gene pairs between pineapple and rice, such as AcoERF15, AcoERF22 and AcoERF32, suggesting that the appearance of these orthologous gene pairs before the divergence of monocotyledonous and dicotyledonous plants. In order to further understand the evolution of ERF gene family, the Ka/Ks ratios of the ERF gene pairs were calculated. The Ka/Ks ratio can represent different selection categories for duplication genes, such as the Ka/Ks >1 indicates positive selection , the Ka/Ks = 1 indicates neutral evolution, and the Ka/Ks <1 indicates negative selection (Li et al., 2019). According to the results, the Ka/Ks ratio of most orthologous ERF gene pairs are less than 1, indicating that they might have experienced strong negative selection pressure during the evolution of pineapple.

Expression patterns of AcoERF genes in different tissues of pineapple

The expression profiles of 74 AcoERF genes in various pineapple tissues at different development stages was performed using RNA-seq expression data from MD2 pineapple plants recently published (Ming et al., 2015). According to the result, a total of 70 AcoERF genes were expressed in different tissues, whereas 4 AcoERF genes (AcoERF22, AcoERF32, AcoERF43, AcoERF72) were not found in the RNA-seq libraries (Fig. 6, Table S6). 22 genes were expressed in all tested tissues, in which 5 genes (AcoERF13, AcoERF16, AcoERF31, AcoERF42, AcoERF65) showed relative high expression levels (value > 10), suggesting these genes might played the indispensable roles in pineapple development. Moreover, 8, 16, 6, 9 genes showed high expression levels at all selected stages in ovule, sepal, stamen and petal, respectively; 16, 13, 12 genes showed high expression levels in leaf, flower and root, respectively; 10 genes showed high expression levels in fruit at all detected stages, indicating that these genes might be involved in the growth of relative tissues. Additionally, AcoERF42 and AcoERF65 showed the highest expression levels in fruit, suggesting that these genes might play positive roles in the growth and development of fruit.

Identification of cis-elements in AcoERF promoters

In order to determine the response mechanisms of AcoERFs, 13 stress-related cis-elements were found in pineapple ERF promoters, such as low temperature response elements (LTR), dehydration responsive elements (DRE) and defense and stress responsive elements (TC-rich) (Table S7). LTRs are involved in low-temperature stress response (Choudhury et al., 2008); DRE and TC-rich repeats play roles in dehydration, low temperature, and salt stress response (Diazdeleon, Klotz & Lagrimini, 1993; Germain et al., 2012). According to the results, all AcoERF genes contain more than two cis-elements in their promoters. The promoter of AcoERF23 only contained two cis-elements, whereas the most ten in the promoter of AcoERF46. In addition, at least one ABRE was present in 85.13% (63 out of 74) AcoERF gene promoters. The analysis of the 13 cis-elements suggested that most of AcoERF genes could response to different stress conditions.

Expression profiles of AcoERF genes under various abiotic stresses

According to previous studies, some ERF genes were involved in abiotic stress response in various plants, such as Arabidopsis (Sakuma et al., 2002), grape (Zhuang et al., 2009), tomato (Sharma et al., 2010), peach (Amygdalus persica L.) (Zhang et al., 2012). However, there was no report about the function of ERF genes in response to abiotic stress in pineapple. To investigate the roles of AcoERF genes in abiotic stress response, 9 AcoERF genes in roots and leaves were selected for functional verification. RT-PCR experiments were performed to analyze the expression patterns in response to different abiotic stress treatments, including cold, drought and salt stress (Figs. 7–9, Tables S8–S10). Overall, we found that the expression of these genes was influenced by the treatments. Among these treatments, some AcoERF genes in different pineapple tissues were significantly induced by the treatments, such as AcoERF36 AcoERF45, and AcoERF67. Interestingly, the expression patterns of AcoERF genes from root and leaf were different under the same treatments, such as AcoERF13, AcoERF16 and AcoERF42 under cold stress, AcoERF7, AcoERF36 and AcoERF60 in response to drought and salt stress, indicated that these genes may have different functions in various tissues of pineapple. Under cold and drought stress, the expression level of most AcoERF genes in leaves was higher than that in roots, while it was the opposite under salt stress, suggesting that AcoERF genes in root were sensitive to salt stress. Besides, most of AcoERF genes had maximal expression before 12 h in drought and salt stress, indicating these genes might play key roles in the early stage of drought and salt stress responses.