Genome-wide analysis of AP2/ERF transcription factors in zoysiagrass, Zoysia japonica

ABSTRACT The AP2/ERF (AP2, APETALA2; ERF, ethylene responsive factor) transcription factors regulate various important processes in plants, including growth, development, metabolism and stress responses. In the zoysiagrass genome (Zoysia japonica), 147 AP2/ERF genes were identified, which were classified into three families based on phylogenetic analysis. By analysis of RNA-sequencing data, the expression profiles of 46 members were assessed, and their regulation patterns in response to cold stress were investigated. Meanwhile, a zoysiagrass genetic regulation network (GRN) was reconstructed based on rice interaction data using the homology search method. Two GRN modules containing AP2/ERF transcription factors were characterized as having important functions in exocytosis and sucrose metabolism processes, which suggested that AP2/ERF transcription factors play critical roles in these biological processes, exerting control over plant growth and stress responses. These results provide useful information for future functional analyses of the AP2/ERF family in zoysiagrass.


Introduction
Plant growth and development are influenced by many biotic and abiotic stresses, such as cold, drought and high salinity. To adapt to the environment, which might be unfavourable for growth, plants employ complex molecular networks, involving many functional genes, such as transcription factors (TFs). The AP2/ERF TFs are one of the largest groups in plants and comprise approximately 60-70 amino acids. They are involved in DNA binding, and are named after their APETALA2 (AP2) domains. Based on the number of AP2 domains, and the presence of other DNA binding domains, AP2/ERF TFs can be divided into the AP2, ERF (ethylene responsive factor), RAV (Related to ABI3/VP) and Soloist families. The TFs of the AP2 family have two AP2 domains and function mainly during flower development. The members of the ERF family have a single AP2 domain, and are divided into two subfamilies based on their binding sites: the ERF and DREB (dehydration responsive element binding) subfamilies. TFs of the ERF subfamily, which bind to GCC-boxes, are involved in hormone signalling pathways, such as the ethylene, jasmonic acid and salicylic acid pathways [1][2][3]. Members of the DREB subfamily bind to CCGAC sites, and are involved in gene expression responses to various abiotic stresses, such as drought, cold and high salinity [4,5]. The RAV family TFs include two domains, an AP2/ERF domain and a B3 domain, and have also been characterized as regulators in plant development and abiotic stress responses [6,7]. The remaining TFs are named as Soloist, with only one member identified in Arabidopsis [8].
To date, AP2/ERF TFs have been identified and investigated in numerous plants, such as Arabidopsis [9], rice [9], poplar [10], soybean [11] and grapevine [12]. Therefore, numerous researches involving the genome-wide analysis of AP2/ERF TFs have contributed to revealing their regulatory functions in plant growth and development, and especially, stress responses; however, they are reported rarely in non-model plants. Zoysiagrass (Zoysia japonica) is a widely used, environmentally friendly, warm-season turf grass species, whose genome has been completed and published recently. As a perennial grass, cold stress is the limiting factor for zoysiagrass distribution in transitional and temperate regions. Therefore, understanding cold tolerance has become a critical issue in zoysiagrass genetic breeding. Wei et al. [13] performed high-throughput sequencing to explore the cold response in zoysiagrass, and identified three AP2/ERF genes that are differentially expressed under cold stress.  However, most members of the AP2/ERF TFs have not been characterized in zoysiagrass. Therefore, in the present study, we performed a comprehensive analysis of AP2/ERF TFs in zoysiagrass, including gene identification, classification and phylogenetic analysis. Their expression profiles under cold stress were also investigated using high-throughput sequencing data. Finally, the genetic regulation network (GRN) of zoysiagrass was reconstructed based on the rice GRN, and AP2/ERF TFs involved in critical regulation relationships were explored using network analysis, which will contribute to revealing the molecular mechanism underlying the involvement of AP2/ERF TFs in growth, development and stress response processes.

Materials and methods
Identification and classification of AP2/ERF genes in zoysiagrass Z. japonica genome sequences were downloaded from the 'Zoysia Genome Database', and the sequence information was listed as described by Tanaka et al. [14] (http://zoysia.kazusa.or.jp). The protein sequences of AP2/ERF genes in rice and Arabidopsis were collected from the Transcription Factor Database (http://planttfdb. cbi.pku.edu.cn/) [15], and these protein sequences were used as query sequences for a BLAST (Basic Local Alignment Search Tool) search against the zoysiagrass genome with the following parameters: expected values 1E-5 and more than 80% coverage. All BLAST hits were retrieved and searched for AP2 domains using the HMMER algorithm, with a cut-off value of 1E-2. The hidden Markov model (HMM) profile of the AP2 domain (PF000847) was downloaded from the Pfam website (https://pfam.sanger.ac.uk) [16]. Genes encoding proteins with AP2 domains were identified as AP2/ERF gene candidates. All the candidates were aligned with the rice and Arabidopsis AP2/ERF genes to classify them into different groups, as described by Nakano et al. [9]. In addition, all the annotation information for the candidate genes was retrieved from the zoysiagrass genome annotation, such as introns and protein length.

Phylogenetic analysis of the zoysiagrass AP2/ERF genes
To explore the phylogenetic relationship of zoysiagrass AP2/ERF genes, multiple alignments of candidate AP2/ ERF genes were performed using ClustalW with default parameters [17]. The results were used to reconstruct phylogenetic trees using MEGA (V4.0) [18]. The phylogenetic trees were generated using complete protein sequences with the following parameters: pair-wise deletion, Poisson correction and 1,000 bootstrap replicates.

Expression analysis of zoysiagrass AP2/ERF genes in response to cold stress
The RNA-seq data of zoysiagrass response to cold stress were downloaded from the NCBI SRA database. The data contained three zoysiagrass samples; the information about the data and the experiment was described in detail in [13]. All the high-throughput sequencing reads were mapped to the zoysiagrass genome using TopHat2 [19], and the gene expressional levels were evaluated using Cufflinks [20], according to the software's instructions. All AP2/ERF gene expression data were retrieved, analysed and clustered using the R program.

Genetic regulation network analysis of zoysiagrass AP2/ERF genes
The rice genome-wide protein-protein interaction network was downloaded from RicePPINet (http://netbio. sjtu.edu.cn/riceppinet), and the genetic GRN was shown as described by Liu et al. [21]. All proteins were BLAST searched against rice proteins, with an e-value cut-off of 1e-05, and hits with highest scores were identified as homologous genes for the zoysiagrass genes. The GRN of zoysiagrass was then reconstructed using the rice GRN data, based on homology between zoysiagrass and rice. The zoysiagrass GRN was analysed using Cytoscape and ClusterViz with the MCODE method and default parameters [22,23]. Sub-networks containing more than three AP2/ERF genes were retrieved and analysed, and the genes of these sub-networks were annotated using Gene Ontology information based on the homology search results. The sub-networks were subjected to gene ontology (GO) enrichment analysis using topGO, and the high enrichment terms were assigned as GRN functions, as described in the software's protocol [24,25].

Results and discussion
Identification, classification and phylogenic analysis of AP2/ERF genes in zoysiagrass Based on homology searches and domain analysis, we identified 147 candidate AP2/ERF genes in the zoysiagrass genome, which were named as ZjERF001 to ZjERF147. These AP2/ERF genes encode putative proteins with lengths of 116 to 654 aa, and have 0 to 9 introns (see Supplementary Table S1). According to the annotations of rice and Arabidopsis AP2/ERF genes, they were assigned to three families: the ERF family (131 members, with a single AP2 domain), the AP2 family (10 members, with two AP2 domains) and the RAV family (six members, with a single AP2 domain and a B3 domain) (Figure 1). The ERF family was divided into two subfamilies based on their homology with members from model plants, including rice and Arabidopsis. The ERF subfamily contained 52 members (ZjERF001-ZjERF052), and was subdivided into six groups (A1-A6). The DREB subfamily with 79 members (ZjERF053-ZjERF131) was subdivided into seven groups (B1-B7), which was consistent with a previous report [9]. The ERF subfamily functions mainly in plant growth, while the DREB subfamily is involved mainly in abiotic stress responses, especially to cold and drought stresses [26][27][28][29][30][31]. There were more members of the DREB subfamily than the ERF subfamily, which could suggest that there are more complex genetic regulations in zoysiagrass, allowing it to grow in unfavourable environments, such as low temperature or water deficiency.

Expression analysis of AP2/ERF genes in response to cold stress in zoysiagrass
To explore the expression profiles of AP2/ERF genes in response to cold stress, we analysed RNA-seq data from a public database. The results showed that 46 ZjERF genes were differentially expressed in response to cold stress and were clustered into four groups (Figure 2). The expression of genes in group A was down-regulated  by cold stress, and comprised 12 members, including three members of the RAV family (ZjERF143, ZjERF145 and ZjERF146). RAV TFs were reported to be negative regulators of plant growth and abiotic stresses, and our results showed their reduced expression levels under cold stress, which supported their negative regulation function for cold tolerance in zoysiagrass [6]. The genes in the other groups (B, C and D) were induced by cold stress; however, they did not show a uniform response to cold stress. The three groups contained 32 AP2/ERF genes. Except ZjERF147, which is a member of the RAV family, most of them were members of ERF or DREB subfamilies, and their expressional profiles showed divergent responses to cold stress. The members of groups B and C, containing 12 members, were highly expressed at 2 h after cold stress, implying that they responded rapidly to cold stress. While group D, comprising 19 members, was highly expressed at three days after cold stress, which indicated a long-term regulatory function in response to cold stress. Previous studies have shown that cold regulation by AP2/ERF genes involves an extensive co-regulation network, including those functioning as first-wave TFs and downstream regulatory TFs [32]. Our findings also suggested that zoysiagrass AP2/ERF genes, especially the members of the ERF and DREB subfamilies, probably participate in a low-temperature regulatory network to determine cold tolerance in zoysiagrass, similar to other plants [32].
Genetic regulation network analysis of zoysiagrass AP2/ERF genes Using the homology search method, we reconstructed the zoysiagrass GRN with 9132 genes (nodes) and 215163 interaction (edges), based on the rice GRN data. To investigate the genetic function of the AP2/ERF genes, we clustered the zoysiagrass GRN, and two modules involving more than three AP2/ERF genes were selected and named as SubModule08 and SubModule69 ( Figure 3). SubModule08 comprised 257 genes and 1212 interactions, and contained six members of the AP2/ERF superfamily (ZjERF009, ZjERF014, ZjERF088, ZjERF102, ZjERF121 and ZjERF127). According to the GO annotation analysis, this GRN module was enriched significantly for the exocytosis process (GO:0006887, p value: 6.6e-08; Supplementary Figure S1), implying that the six AP2/ERF genes played important roles in exocytosis. In melon, Corbacho et al. [33] investigated ERF TFs with critical functions in controlling fruit development via regulation of exocytosis, which was also identified in tomato. The findings concerning SubModule08 were consistent with previous studies, implying that zoysiagrass AP2/ERF genes would also control plant growth by regulating exocytosis [33]. SubModule69 contained 64 genes and 102 interactions, with four AP2/ERF genes (ZjERF007, ZjERF023, ZjERF095 and ZjERF107). This GRN module was enriched for the sucrose metabolic process (GO:0005985, p value: 3.0e-07; Supplementary Figure  S2). Sucrose is an important osmotic adjustment compound, which has been characterized widely and has critical roles in abiotic stress in numerous plants, such as Arabidopsis [34], rice [35] and poplar [36]. These studies showed that AP2/ERF TFs regulate sucrose metabolism and modify the response to abiotic stress. In the present study, four members of the zoysiagrass AP2/ ERF family were characterized as being involved in sucrose metabolism, implying that a similar regulatory pathway is present in zoysiagrass to confer abiotic stress tolerance.

Conclusions
In the present study, 147 members of the AP2/ERF superfamily were identified and characterized in zoysiagrass. The expression profiles of 46 members were assessed from RNA-seq data in response to cold stress, which confirmed their regulatory functions in the cold stress response. Meanwhile, a zoysiagrass GRN was reconstructed using rice interaction data, and two GRN modules involving AP2/ERF TFs were analysed in detail. The results indicated that zoysiagrass AP2/ERF TFs play important roles in plant growth and stress response by regulating exocytosis and sucrose metabolism. Our findings will provide valuable information to determine the functions and molecular mechanisms of AP2/ERF TFs in zoysiagrass.

Disclosure statement
No potential conflict of interest was reported by the authors. Note: The yellow nodes represent ZjERF genes; the pink nodes represent other genes and the green edges represent regulation relationships between genes in zoysiagrass.