Previously, we established a comprehensive transcriptome of S. caninervis transcripts during a standardized dehydration-rehydration process (Gao et al., 2014). Ninety-two thousand, two hundred and forty S. caninervis unigenes were obtained and served as the source for the AP2/ERF genes identification presented in this study. Arabidopsis and Physcomitrella patens AP2/ERF family genes were downloaded from the plant TF database (PlantTFDB v3.0)¹ (Jin et al., 2014). The HMM profiles were downloaded from Pfam database 27.0² (Punta et al., 2012).

One hundred and seventy-six Arabidopsis AP2 predicted amino acid sequences and one hundred and seventy-one P. patens AP2 predicted amino acid sequences were used as queries to search against the S. caninervis transcriptome database using BLASTP programmer, and an E value of 1e-3 was adopted. Next, the HMM profile PF00847 (AP2 domain) and PF02362 (B3 domain) were queried using hmm search command included in the HMMER (v3.0) software with an E value cutoff at 1e-3. All the candidate ScAP2/ERF genes identified through these two methods were confirmed with Conserved Domain Database (CDD³; Marchler-Bauer et al., 2015) and SMART⁴ (Letunic et al., 2015) search to ensure the presence of AP2 domain. An AP2 domain length of approximately 60 amino acids was considered to be a full-length AP2 domain (Duan et al., 2013). Only sequences longer than 200 nucleotides were retained for further analysis. Sequences which shared >95% matches were considered redundant. The detailed flowchart of AP2/ERF gene detection in S. caninervis is depicted in Supplementary Figure S1.

Open reading frames (ORFs) were predicted with the ORF Finder at NCBI⁵. Motif detection was performed with the online tool MEME program⁶ (Bailey et al., 2009) using the parameters: any number of repetitions per sequence, motif width ranges of 6–50 amino acids, and 11 as the maximum number of motif. The motif alignment was conducted using MAST⁷ and Patmatch programs⁸ (Bailey et al., 2009). Physical and chemical characterization of genes with deduced amino acid sequence were analyzed using ProtParam⁹ and protein subcellular location(s) were predicted using PSORT¹⁰. Multiple sequence alignment was performed with ClustalW in conjunction with MEGA 5.1 (Tamura et al., 2011), phylogenetic trees were constructed by the neighbor-joining (NJ) method using MEGA 5.1 (Poisson correction and pairwise deletion). Support for nodes on the estimated phylogeny was tested with 1000 bootstrap replicates.

Statistical analyses were performed using SPSS software (standard version 11.5 released for Windows, SPSS Inc., Chicago, IL, USA). All data were analyzed using a one-way analysis of variance (ANOVA) at the 95% confidence level. The significant difference was compared using LSD multiple comparison test. The data shown are the mean values ± SE of three replicates, and the significance level relative to controls is ^∗P < 0.05, ^∗∗P < 0.01. Figures were generated using Sigmaplot10.0. Adobe Illustrator CS5 and Adobe Photoshop CS3 were used for image processing.

Eighty AP2/ERF unigenes were identified from the S. caninervis transcriptome using HMMs and BLASTP (see Supplementary Figure S1 for the detailed screening strategy). The unigenes ranged from 202 to 2514 bp in length, and the corresponding predicted amino acid sequences ranged from 39 to 659 aa. AP2 domain prediction of the ScAP2/ERF genes demonstrated that 68 of the 80 unigenes (85%) contained a full-length AP2 domain (ca 60 aa), five unigenes contained two AP2 domains, and a single unigene contained both a single AP2 domain and a B3 domain. Seventeen of the 80 unigenes (21%) contained an intact ORF. Based on these results, the 74 unigenes encoding a single AP2 domain which were preliminarily classified as members of the ERF protein family (some special AP2 family protein members and Soloist proteins also contain a single AP2 domain). The five AP2 double-domain unigenes were classified as AP2 protein family members and the AP2/B3 unigene was classified as a RAV protein family member. Group classification of the ScAP2/ERF genes was accomplished by constructing a phylogenetic gene tree which aligned the predicted AP2 domains using the NJ method. Nine unigenes were excluded from further analysis due to truncated AP2 domains, finally 71 predicted amino acid sequences were used to construct the phylogenetic tree. The un-rooted gene tree demonstrated that the 71 S. caninervis AP2/ERF predicted amino acid sequences can be divided into five main groups: ERF I, ERF II, ERF III, AP2, and RAV (Supplementary Figure S2).

Sixty-five unigenes containing a single AP2 domain grouped into three clades: ERF I (two sequences), ERF II (41 sequences), and ERF III (22 sequences). Five unigenes which contained two AP2 domains grouped into the same clade (AP2), and the AP2/B3 unigene formed a unique clade (RAV). Additional phylogenetic analysis employing 71 ScAP2/ERF predicted amino acid sequences and 176 previously annotated Arabidopsis AP2/ERF predicted amino acid sequences supported the classification into five groups (Figure 1). The number of ScAP2/ERF subfamily genes were compared with the model plants Arabidopsis and P. patens. The ERF subfamily was the largest in each species ranging from 77 to 90% (Table 1), while RAV and Soloist subfamily was the smallest ranging from 1 to 4%.

Within the plant AP2/ERF family, the ERF subfamily is dominant and has been extensively documented to play important roles in plant responses to biotic and abiotic stresses (Zhuang et al., 2008; Huang et al., 2016). Phylogenetic analysis using 63 ScERFs and 136 AtERFs (both of which excluded soloists) identified members of eight out of 12 subfamilies based upon the classification of Sakuma et al. (2002): A2, A3, A4, A6, B1, B3, B4, and B6. ScERF predicted amino acid sequences contained no members of the A1, A4, B2, and B5 subfamilies (Figure 2 and Supplementary Figure S3). Accordingly, based upon Nakano et al. (2006)’s classification method, ScERF genes contained members of group I to group Xb-L, with the exception of group II, III and group VI. Most plants have a single A-3 type DREB (Sakuma et al., 2002; Zhang et al., 2008) while some species (such as grape) lack A-3 type of DREBs (Licausi et al., 2010). Unigene26481 clustered together with both A-2 and A-3 type of AtDREBs (Figure 2) which prevented a definitive classification. Additional analysis using 57 AtDREBs generated a tree that showed unigene26481 was clearly grouped together with the A-3 type DREB in Arabidopsis (Supplementary Figure S3). S. caninervis has a single A-3 DREB protein. Interestingly, more than half of the ERF sequences are members of the A-5 subgroup DREBs in S. caninervis while other subgroups contain between two and six sequences (Figure 2 and Supplementary Figure S3).

To analyze the conservation of the AP2 domains, ScERF deduced amino acid sequences were aligned with 22 AtERF sequences representing each ERF subfamily/group. Sequence alignment demonstrated that ScERF sequences shared significant amino acid similarity with AtERFs and that the AP2 domains of ScERFs also contained three β-sheets and one α-helix (Figure 3). The amino acid residues 4 G, 16 E, 27 W, 28 L, 29 G were completely conserved in all ScERF and AtERF genes (labeled with asterisk). β-sheet 1 (β-1) contains the conserved GVR element (4G, 5V, 6R) with several sequences substituting 5 I for 5V. β-sheet 2 (β-2) contains the conserved EIR element (16 E, 17 I, and 18 R). We found that in both Arabidopsis and S. caninervis, the amino acid residues EVR were observed only in A-1 group DREB genes, and the residues ERK were observed only in the B-6 protein subfamily (Figure 4). These conserved amino acid patterns may be helpful for the classification of ERF genes in other species. In addition, the unigene49249 in S. caninervis had a particular “EMR” element in this position which was not found in Arabidopsis ERF genes. β-sheet 3 (β-3) contained the conserved WLG element (27 W, 28 L, 29 G) which was highly conserved between Arabidopsis and S. caninervis, and the α-helix contained a consensus sequence AAxAxD which was conserved between Arabidopsis and S. caninervis (Figure 3 and Supplementary Figure S4).

Seventeen out of 80 ScAP2/ERF genes (2 ScAP2, 2 ScSoloists, and 13 ScERFs) were predicted to encode a single, intact ORF. The predicted polypeptide was characterized and the putative subcellular localization was analyzed (Table 2). The 17 ScAP2/ERF predicted polypeptides ranged from 172-to-371aa residues. The predicated molecular mass of the deduced polypeptides ranged from 19.13 to 42.27 KDa, and the theoretical pI values ranged from 4.62 to 10.0. More than half of ScAP2/ERF deduced polypeptides were predicted to localize in the nucleus (including two ScAP2). Four deduced polypeptides localized to the cytoplasm and one each localized to the microbody and plasma membrane. The 17 ScAP2/ERF predicted ORFs were submitted to MEME analysis to investigate the motif composition of the deduced polypeptides. In total, 11 motifs were found and their amino acids sequences were showed in Figure 4B. Five out of eleven motifs (motif 1, 2, 3, 7, 8) represented the AP2 domain. Motif 2 and 7 correspond to the AP2 domain β2-and β1-sheet of AP2 domain, respectively, and Motif 1 contained the β3-sheet and α-helix. Motif 3 was similar to Motif 8 which corresponded to the C-terminal of the AP2 domain (Figure 4). ScERF genes with intact ORF should be the priority candidates for gene cloning and function analysis. Therefore, we constructed the gene tree using 15 ScERF genes (including two Soloists) with intact ORFs and 139 Arabidopsis ERF genes to establish a detailed classification. The tree revealed that 15 ScERFs covered two ScSoloists, 12 A-5 type of DREBs, and one B-4 type of ERF gene (Figure 5). Three genes (unigene 40134, 48763, and 86206) belonged to A-5a subgroup, two genes (unigene 9249 and 7018) belonged to A-5b subgroup, and unigene 74743 which had the longest ORF belonged to A-5c subgroup. Additionally, six other A-5 type ScERF genes cannot be classified into a subgroup a, b nor c, which separately grouped in a S. caninervis-specific clade. One gene (unigene83605) was classified into the B-4 group.

Based on the phylogenetic tree, eight ScDREBs representative of different group (such as A-5a, A-5b, A-5c) were selected to quantify the changes in transcript abundance during the dehydration-rehydration process using RT-qPCR analysis. RT-qPCR results showed that compared to full-hydration (labeled as 0 h), all eight transcripts accumulated both during the fast dry (FD) and rehydration (FDR) process (Figure 6). Under FD treatment, transcript abundance peaked at 6 h and then decreased, with the exception of unigene 7018 which peaked at 12 h (Figure 6). During FDR process, transcript abundance increased at 0.5 h and subsequently decreased, with the exception of unigenes 80030 and 9249. The majority of the tested genes have higher transcript abundance during FD compared with FDR. For example the expression level of unigene 9249 increased 12-fold at FD 6 h, but was less than 4-fold increased during rehydration, suggesting that these genes may play a dominant role in fast drying process. While for unigene 40134, which has higher expression level during rehydration process, the gene expression level increased up to 12-fold at FDR 0.5 h compared to the control, suggesting unigene 40134 may play important roles at the early stage of rehydration.

Phylogenetic tree analysis showed the suitability of using the classic AtAP2/ERF-based classification method with moss derived sequence data and provides evidence that the phylogenetic topology of the AP2/ERF family had already been established before the divergence of vascular plants (Mizoi et al., 2012).

The AP2 domain of AP2/ERF genes was reported to be highly conserved among different plant species, although the sequence similarity outside the AP2 domain was very low (Kizis et al., 2001). Similar to other plants, our results showed that the amino acid composition of the AP2 domain was very conservative between Arabidopsis and S. caninervis, and the 4 G, 16 E, 27 W, 28 L, 29 G amino acids were completely conserved in all ScERFs and AtERFs. Other amino acid variation can further refine the classification of genes. For example, the motifs “HLG” and “WLG” in the β3-sheet can distinguish a Soloist gene from an ERF gene. Interestingly, we found that in β2-sheet of AP2 domain of both S. caninervis and Arabidopsis, the vast majority of genes shared an “EIR” element, while “EVR” existed only in the A-1 group of DREB genes and “ERK” was specific to the B-6 subfamily genes in this position.

ScAP2/ERFs showed commonalities as well as unique characteristics as compared with Arabidopsis. Most of ScERF genes can be classified as described above. A small number of ScERF genes form a unique clade (A-5 DREB) but cannot be classified into a subgroup. The AP2/ERF superfamily in moss encompassed all of the Arabidopsis, and the ratio for each subfamily was constant between S. caninervis and P. patens. However, compared with Arabidopsis, the percentage of ERF subfamily members in the moss species were greatly expanded, and the gene members in the AP2 and RAV subfamilies were reduced. This adjustment of subfamily members may connected with their proposed functions in both development and stress response. The AP2 subfamily is associated with plant development and ERFs play important roles in abiotic stress tolerance regulation. Abundant ERF genes in moss might play important roles in abiotic stress tolerance, which is an adaptive evolutionary adaption for plant species growing in adverse environments (Licausi et al., 2013). Genes from different group may have specific motifs. Conserved motifs outside the AP2 domain region may function as a repression or activation domain, such as ERF-associated amphiphilic repression (EAR: DLNxxP) motif (Ohta et al., 2001; Kagale and Rozwadowski, 2011), LWSY motif (Dubouzet et al., 2003; Chen et al., 2008), and DELL motif (Tiwari et al., 2012). In ScERFs, unigene 74743 had the longest ORF and the AP2 domain was located in the end of C-terminus and contained motif 9 (SXSZPISEDSY) which was repeated five times before the AP2 domain (Figure 4). Bioinformatic analysis has failed to identify ‘motif 9’ in any additional plant protein. Motif 9 is a novel sequence and the function of this motif is unclear. A priority of our future work is the functional analysis of this novel sequence.

Dehydration-responsive element-binding proteins genes can be divided into the A1–A6 subgroups according to the sequence similarity of AP2 domain. In Arabidopsis, A-1 (DREB1) and A-2 (DREB2) type DREB proteins have been intensively studied and were reported to play important roles in plant response to cold and osmotic stresses (Sakuma et al., 2002). DREB A-1 and DREB A-2 homologous have been identified and functional analyzed in a variety of plants (Sakuma et al., 2006; Qin et al., 2007; Cong et al., 2008; Yang et al., 2011; Mizoi et al., 2013). Despite representing a large proportion of the DREB family, A-5 type of DREB proteins are poorly studied, and the function and stress response mechanisms of A-5 subgroup proteins is still unclear. In Arabidopsis, stress-associated DREB A-5 proteins were the largest subgroup (with 16 gene members; Sakuma et al., 2002), only the RAP2.1 gene was studied in detail (Dong and Liu, 2010). Additional A-5 DREBs were reported in cotton (Huang and Liu, 2006), soybean (Chen et al., 2007, 2009), potato (Bouaziz et al., 2012), Malus sieversii (Zhao et al., 2012), and Halimodendron halodendron (Ma et al., 2015). Based upon the literature, A-5 DREB proteins showed a diversity in gene expression (i.e., stress response) and putative function (i.e., stress tolerance). The A-5 type of DREBs response to at least one kind of abiotic stress, such as soybean GmDREB3 responds to a single abiotic stress (i.e., low temperature), while GmDREB2 and PpDBF1 in P. patens can respond to drought, salt, cold, and ABA treatments (Chen et al., 2007, 2009; Liu et al., 2007). Moreover, manyA-5 type of DREB proteins have been shown to enhance stress tolerance such as GmDREB2, GmDREB3, and MsDREBA5 (Chen et al., 2007, 2009; Zhao et al., 2012). Interestingly, RAP2.1 in Arabidopsis was reported to negatively regulate drought and cold tolerance in transgenic Arabidopsis (Dong and Liu, 2010). P. patens has a single A-5 type DREB (PpDBF1) isolated using FDD technology (Liu et al., 2007). In this study, the vast majority of DREB genes found in S. caninervis transcriptome were A-5 type DREBs, and the majority respond to dehydration and/or rehydration treatment. Unlike Arabidopsis and other plants, we postulate that A-5 type DREB proteins in moss species are the dominant regulator for stress response rather than A-1 (DREB1) and A-2(DREB2) proteins. Hence, to better understand the functional mechanism of plant response to stresses, more members of A-5 subgroup DREB proteins needed to be identified and studied.