Molecular cloning and expression analysis of three ThERFs involved in the response to waterlogging stress of Taxodium ‘Zhongshanshan406’, and subcellular localization of the gene products

As a subfamily of the APETALA 2/ethylene response element binding protein (AP2/EREBP) transcription factor superfamily, the ethylene response factor (ERF) is widely involved in the regulation of growth and response to various abiotic stresses in plants, and has been shown to be the main transcription factor regulating transcription of the genes related to hypoxia and waterlogging stress. In this study, three ThERF genes, with significant differences in expression profile in response to flooding stress, were identified from the transcriptomics data acquired from Taxodium hybrid ‘Zhongshanshan 406’ (T. mucronatum Tenore × T. distichum (L.) Rich) under waterlogging stress: ThERF15, ThERF39 and ThRAP2.3 (GenBank ID: KY463467, KY463468 and KY463470, respectively).The full-length cDNA of each of the three ERFs was obtained using the RACE (rapid amplification cDNA ends) method, and all three were intron-free. Multiple protein sequence alignments indicated that ThERF15, ThERF39 and ThRAP2.3 proteins all had only one AP2-ERF domain and belonged to the ERF subfamily. A transient gene expression assay demonstrated that ThERF15, ThERF39 and ThRAP2.3 were all localized to the nucleus. Real-time quantitative PCR (qPCR) revealed that the expression of ThERF15, ThERF39 and ThRAP2.3 exhibited significant differences, compared with the control, in response to two levels of flooding treatment (half-flooding or total-submergence) of ‘Zhongshanshan 406’. Quantification of ethylene concentration revealed that ethylene was more relevant to the level of expression than the period of flooding treatment. Based on the experimental results above, ThERF15, ThERF39 and ThRAP2.3 were identified as being related to the regulation of downstream flooding- responsive gene expression in ‘Zhongshanshan 406’. ThRAP2.3 is most likely to be a key downstream-response ERF gene to respond to the output of the ethylene signal generated by flooding stress.


INTRODUCTION
In China, there is great potential for the development of afforestation in water-logged areas, due to the abundant resources of shallows and wetlands. Further research on tree species with high waterlogging tolerance would be beneficial to improve the utilization of land resources in coastal areas and wetlands in China. Afforestation is also an important measure to solve the problems of flood disaster, biodiversity decline and wetland restoration along the coastal and riverside areas of China (Cao et al., 2010).
To date, compared with herbaceous plants, fewer xylophyta have been reported to exhibit long-term flooding tolerance (Wang et al., 2016;Yu et al., 2012). The genus Taxodium is well known for its waterlogging tolerance and landscape values (Conner, Mihalia & Wolfe, 2002;Souther & Shaffer, 2000). Interspecific hybridization in Taxodium Rich can combine the best characteristics of superior parents (Zhou et al., 2010). Taxodium hybrid 'Zhongshanshan406', which is an elite clone selected from the hybrid Taxodium mucronatum ♀×Taxodium distichum ♂, developed by Institute of Botany, Jiangsu Province and Chinese Academy of Sciences, showed great improvement in flooding tolerance, and has been widely planted in China (including Dianchi Lake in Yunnan Province, Chaohu City in Anhui Province, Three-Gorge reservoir in Chongqing Province) as a wetland species (Han & Shan, 2012;Ma et al., 2011). Related physiological and transcriptomics studies have demonstrated that 'Zhongshanshan406' is an ideal model plant with which to research the waterlogging characteristics of woody plants (Hua et al., 2017;Qi et al., 2014).
To cope with abiotic stresses, including waterlogging stress, transcription factors (TFs) play a central role by regulating expression of downstream stress-responsive genes via sequence-specific binding to cis-acting elements in the promoters of target genes (Hussain, Kayani & Amjad, 2011;Mizoi, Shinozaki & Yamaguchi-Shinozaki, 2012). As a large superfamily of plant-specific transcription factors, the APETALA2/ethyleneresponsive binding factor (AP2/EREBP) family is involved in a myriad of regulatory processes, such as growth and development, metabolic regulation, and response to a variety of biological and abiotic stresses (Upadhyay et al., 2013). The AP2/EREBP family includes four major subfamilies: the AP2 (APETALA2), DREB (dehydration-responsive element-binding protein), RAV (related to ABI3/VP1) and ERF (ethylene-response factor) subfamilies (Zhuang et al., 2009). Of these, many ERF subfamily members have been isolated and characterized with respect to their major effects, related to responses to plant hypoxia and flooding stress. These ERFs function by regulating gene expression of downstream stress-tolerance genes via the cis-acting ethylene-responsive element (ERE), known as the GCC box (AGCCGCC;Xu et al., 2007). All the ERF subfamily members have a conserved 58-59 amino acid AP2/ERF domain, which makes them capable of performing the above adjustment functions (Fujimoto et al., 2000;Mizoi, Shinozaki & Yamaguchi-Shinozaki, 2012).
Ethylene is a key upstream regulatory component of the ethylene signal transduction pathway, its effects being mediated by signal transduction components, including ERF transcription factor (TF) families (Yin et al., 2012). Ethylene plays a critical role in myriad developmental programs and fitness responses to pathogens and abiotic stress factors (Alonso & Stepanova, 2004;Guo & Ecker, 2004). Many studies have shown that a series of ethylene receptors (ETR1, ETR2, ERS and EIN4), located in the endoplasmic reticulum membrane, can interact with the ethylene molecule to exert inhibitory effects on the kinase activity of CTR1 (constitutive triple response 1), to release the key positive regulator EIN2(ethylene insensitive 2) of the ethylene signaling pathway, and then send the signals toEIN3/EIL1transcriptional regulatory factors and downstream ERF -responsive genes, thus completing the output of the ethylene signal (Bisson & Groth, 2011;Ju et al., 2012;Li, Ma & Guo, 2013).
Despite numerous reports on the ERF subfamilies in herbaceous model plants (Joo et al., 2013;Licausi et al., 2010a;Licausi et al., 2010b), there has been little progress on the isolation and characterization of the ERF genes operating in response to flooding and hypoxia stresses in woody plants, with no reports from Taxodium hybrid 'Zhongshanshan'. Our objectives of this work were to explore the molecular mechanisms of the ThERF s to tolerate waterlogging stress and it's vital for us because it's a research gap in 'Zhongshanshan'. In this study, we identified and characterized three genes encoding ERF proteins (ThERF15, ThERF39 and ThRAP2.3) from 'Zhongshanshan 406', and explored their expression profiles, in comparison with the ethylene-accumulation mechanism in waterlogging stress, hoping to enrich the genetic resources of the ERF TFs in 'Zhongshanshan 406', and to make a preliminary exploration of the molecular mechanisms involved in the accumulation of ethylene at the protein level.

Plant materials
One-year-old cuttings of 'Zhongshanshan 406' plants were obtained from Institute of Botany, Jiangsu Province and Chinese Academy of Science, Nanjing, China. Prior to the experiment, about 100 healthy 'Zhongshanshan 406' plants were carefully transplanted into plastic pots containing 3:1:1 (v/v/v) clay, vermiculite and perlite in July 2016. Each pot was 20 cm in diameter and 25 cm in height. All plants were irrigated fully every two days. A pot tray was placed under each pot to avoid water loss and soil erosion. The plants were allowed to acclimate to the local conditions for 60 days before the treatments were imposed.
In early September, 90 plants were selected on the basis of uniformity of size and development. One of three treatments was applied to each batch of thirty plants: nonflooding (control, CK), half-flooding (HF) and total-submergence (TS).The plants with treatments HF and TS were put into outdoor concrete tanks (2.8 m × 1.7 m each) filled with tap water to different depths. CK plants were placed near the tanks and watered normally. HF plants were flooded to a water level 5 cm above the soil surface. TS plants were completely submerged. On the first, third, fifth, seventh, and ninth days of the treatments, leaf tissue and root tissue (for treatments CK, HF and TS) were sampled for analysis at each of these five time points, the tissue being immediately snap-frozen in liquid nitrogen and stored at −80 • C until extraction of RNA and DNA. Three independent biological replicates were sampled at each time point for each treatment.

Determination of ethylene content
Tissues and time points for sampling were the same as above. For each time point, the leaf tissue and the root tissue were sampled from the three treatments and five biological replicates were prepared for each sample. The tissue in each replicate weighed 0.2 g, and ethylene was extracted from the tissue and assayed using the Plant ETH ELISA Kit (HCB, Vancouver, Canada). The concentration of ethylene in each replicate was determined by comparing the absorbance of the samples to a standard curve.

RNA extraction and cDNA synthesis
The frozen plant tissues for RNA extraction were ground into a fine powder using a mortar and pestle. Total RNA was extarcted from the leaves or roots from the CK, HF and TS treatments at the five sampling time points, using the RNeasy R Plant Mini kit (Qiagen, Dusseldorf, Germany), and was reverse transcribed with iScript TM cDNA Synthesis Kit (BIO-RAD, Hercules, CA, USA), using 1 mg RNA as the template.

Molecular cloning
On the basis of the sequence fragments obtained from the transcriptomics data, which was obtained from 'Zhongshanshan 406' plants under waterlogging stress, nested primers were designed using the Oligo software (Version 6.0) to amplify the full-length sequences with the 3 -Full RACE Core Set Kit (TaKaRa, Otsu, Japan) and SMATer RACE 5 /3 Kit (Clontech, Palo Alto, CA, USA), according to the manufacturer's instructions. The amplified fragments were separated on 1% agarose gels and purified by QIAquick Gel Extraction Kit (QIAGEN, Dusseldorf, Germany), linked into pMD19-T vectors (TaKaRa, Otsu, Japan) and finally transformed into Escherichia coli strain TOPIO. The recombinant plasmids were checked by PCR, and the positive colonies were sequenced. The overlapping sequences were assembled by BioEdit software (Version 2.6) and the full-length cDNA sequences of the three ERFs were obtained. The predicted open reading frames (ORFs) were subsequently amplified by PCR, and were verified by sequencing. Genomic DNA was extracted using a Plant Genomic DNA Extraction Kit (BioTeke, Beijing, China). Genomic DNA sequences of the above genes were amplified with the RNase-treated DNA, and were verified by sequencing. The sequences of the primers are listed in Table 1.

Bioinformatics and statistical analyses
The online BLAST software (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to analyze the DNA and protein sequences of the three ERFs. ORFs were predicted by DNAMAN software (Version 8.0). The physicochemical properties and amino acid composition of the proteins were predicted and calculated using Expasy Protparam (http: //web.expasy.org/protparam/). Analyses of the signal peptide cleavage site, protein domain search and transmembrane structures of the genes were carried out with SignalP online tools (http://www.cbs.dtu.dk/services/SignalP/), PROSITE (http://prosite.expasy.org/) and TMHMM (http://www.cbs.dtu.dk/services/TMHMM/), respectively. GOR IV secondary structure prediction method (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page= npsa_gor4.html) was used to predict secondary structures of the deduced amino acid sequences. Alignment of the deduced protein sequences of ThERF15, ThERF39  (Kumar, Stecher & Tamura, 2016). Data of ethylene concentration approximated to normality, and was analyzed with one-way analysis of variance (ANOVA) followed by multiple comparisons using Duncan's multiple range test at P = 0.05, using SPSS 19.0 software (SPSS Inc., Chicago, IL, USA).

Real-time qPCR analyses
To quantify the expression of ThERF15, ThERF39 and ThRAP2.3, qPCR was performed on an Analytik Jena qTOWER2.2 PCR System (Biometra, Gottingen, Germany). The primers for the three ERFs (Table 1) for qPCR were designed, based on the sequence of their cDNA. Adenine phosphor ribosyl transferase (APRT, GenBank accession No. KX431853) was used as a reference gene, which was amplified using the primers APRT-F and APRT-R (Table 1).
Each sample was carried out in triplicate. The results are displayed in the form of relative expression values 2 − Ct , where Ct represents the Ct value of the gene minus that of the reference gene (Bustin et al., 2009;Livak & Schmittgen, 2001).

Ethylene concentration
Compared with the HF and TS treatment groups, the concentration of ethylene in the CK group fluctuated within a relatively narrow range during the whole treatment process. In both roots and leaves, ethylene production was higher in the TS treatments, compared to the CK, whereas ethylene production in HF and CK were generally not significantly different in either roots or leaves. In both leaves and roots, the lowest ethylene concentrations occurred in the HF treatment on day 9, whereas such a decline was not observed in TS (Fig. 1).
To explore the evolutionary relationship between ERF proteins from different species, the NJ phylogenetic tree was constructed using ThERF15, ThERF39 and ThRAP2.3 sequences together with the other 81 ERF proteins reported from different plants in the  GenBank database, including 65 AtERF proteins (from Arabidopsis thaliana) which have been identified and described in detail (Nakano & Suzuki, 2006) (Fig. 3). The un-rooted NJ tree showed that the 84 proteins were clustered into eleven distinct groups (Fig. 3).
ThRAP2.3 was positioned in the ERF-VII group, while ThERF15, which was most closely related to TwERF15, was positioned in the ERF-IX group, and ThERF39 was located in the ERF-VII-L (VII-Like) group, due to its relative similarity to the ERF-VII group.

Expression profiles of the ERF genes
To analyze the expression profiles of ThERF15, ThERF39 and ThRAP2.3 in 'Zhongshanshan 406' under two levels of flooding stress, we carried out qPCR to measure transcript levels at five waterlogging time points, namely 1 d, 3 d, 5 d, 7 d and 9 d, and in different tissues, namely leaves and roots. CK values, being references for two treatments on the same stage and tissue, were set to 1. The results showed that these three ERF genes were expressed, at different levels, at each of the five waterlogging time points under HF and TS treatments in roots and leaves of 'Zhongshanshan 406'. Except for ThERF15 in leaves, the expression levels of all three ERFs exhibited a decline (at around day 3 or 5) during the early treatment period, before exhibiting a rise in the late stage of the TS treatment, in both roots and leaves, generally with a maximum expression on the 9th day. Unlike the similar expression patterns in the two tissues under TS treatment, expression patterns in the roots and leaves were significantly different in the HF treatment. ThERF39 and ThRAP2.3 both exhibited lower expression levels, compared with CK treatment apart from expression of ThERF39 on the 9th day, which showed a marked increase. In contrast with the two genes described above, the expression pattern of ThERF15 had a greater difference in the late period of HF treatment. Significant increases in ERF gene expression generally occurred only on days 7 and 9 (Fig. 4).

DISCUSSION
In the present study, three ERFs from Taxodium hybrid 'Zhongshanshan 406', with differential expression under waterlogging stress, were isolated and characterized. The deduced proteins of these three ERFs contained a typical AP2/ERF domain and hence are considered to belong to the plant-specific ERF subfamily (Li et al., 2017). Multiple sequence alignments revealed that ThERF15, ThERF39 and ThRAP2.3 exhibited considerable protein sequence similarity with other ERFs from Arabidopsis and Populus, including the AP2-ERF domain, the YRG element, the WLG motif, and the nuclear localization signal (NLS) element. Associated with their role as transcription factors (TFs), all three ThERFs have basic amino acid regions that possibly function as NLSs to target the proteins to the nucleus. The nuclear localization was further verified by transient expression of ThERF15, ThERF39 and ThRAP2.3 in Populus protoplasts.
In the meantime, the phylogenetic analysis of plant ERF proteins indicated that ThERF15, ThERF39 and ThRAP2.3 were classified into the ERF-IX, ERF-VII-L and ERF-VII groups, respectively. Studies have shown that ERF proteins act as regulators of plant development and response to biotic and abiotic stresses, including drought, salt and flooding (Quan et al., 2010). Members of the ERF-VII group were confirmed to respond to plant hypoxia and waterlogging stresses and to regulate the expression of related downstream genes in Arabidopsis and rice (Hattori et al., 2009;Licausi et al., 2010;Licausi et al., 2010b;Xu et al., 2006), whereas ThRAP2.3, which is a member of the ERF-VII group in 'Zhongshanshan 406', is extremely likely to be involved in regulation of the waterlogging stress response. However, there are few reports on the flooding-resistant function of members of the ERF-IX and ERF-VII-L groups, and related research needs to be carried out in detail. Transcript levels of ThERF15, ThERF39 and ThRAP2.3 were measured by real-time PCR at five waterlogging time points under the HF and TS treatments. Expression patterns of the three ERF genes showed that all three exhibited a sharp increase in expression levels in the latter stages of TS treatment, reaching peak expression on the 9th day, except for ThERF15 in leaves. Contrary to the results from the TS treatment, the three ERFs showed varying expression profiles under the HF treatment. ThERF39 and ThRAP2.3 generally had lower expression levels compared with the CK treatment except for ThERF39 on the 9th day. The expression pattern of ThERF15 showed greater differences in the later stages of HF treatment. Despite the different expression patterns of the three ERFs under the same treatment, all of them exhibited significant increases in expression during the time-course of the two flooding treatments. This indicates that the expression of ThERF15, ThERF39 and ThRAP2.3 are clearly influenced by flooding stress, and that the three ERF genes may very possibly be involved in response to flooding stress, especially to extreme flooding stress. The gene-specific response mechanisms and the fine regulation of expression of ThERF15, ThERF39 and ThRAP2.3 deserve further study. In general, the expression of ERF genes is in relation to the molecular response to ethylene. Many ERF TFs are indeed ethylene responsive, despite the observation that the AP2/ERF TFs to which they belong are regulated by numerous physical-chemical stimuli (Licausi et al., 2010;Licausi et al., 2010b). Ethylene has been verified to play a pivotal role in plant responses to biotic and abiotic stresses, including flooding stress (Bleecker & Kende, 2000;Voesenek et al., 2016;Zhang et al., 2016), and these responses are possibly mediated through the ERF gene family, which serve as regulatory elements (Yin et al., 2012). In this study, the concentration of ethylene was determined in the same tissues and treatments as used for the gene expression studies, to explore the regulatory mechanism between plant ethylene accumulation and expression of the three ERFs in 'Zhongshanshan 406'.
Under the HF treatment, the ethylene concentration differed little from that in the CK treatment, although it exhibited an obvious decline during the latter stages of the time-course. During the time-course of the TS treatment, the ethylene concentration in plants under the flooding treatment was always higher than that in plants under the CK treatment. By comparing the expression levels of the three ERF s and the trends in ethylene accumulation during flooding stress, it could be seen that ethylene accumulation in roots and leaves showed a clear increase in the early stages of the TS treatment, while the transcription level of ThRAP2.3 at this point was generally low, increasing, in both leaves and roots, in the later stages of the TS treatment. On the other hand, in comparison with the CK treatment, the accumulation of both ethylene and ThRAP2.3 mRNA showed a similar trend under the HF treatment, both parameters exhibiting a decline in the later stages of the HF treatment. These findings indicated thatThRAP2.3 is most likely to be one of the key downstream-response ERF genes to respond to the output of the ethylene signal generated by flooding stress. No direct connection between ethylene content and the transcription level of either ThERF15 or ThERF39 could be found. There is a high probability, therefore, that expression of ThERF15 and ThERF39 was coordinated by a combination of ethylene and other phytohormones, such as jasmonate (Lorenzo et al., 2003;Sherif et al., 2012). The mechanism of response of expression of ThERF15 and ThERF39 to flooding stress, therefore, needs further research.

CONCLUSIONS
Taken together, these results strongly indicate that ThERF15, ThERF39 and ThRAP2.3 play an essential part in tolerance to flooding stress in 'Zhongshanshan 406', especially ThRAP2.3. The combined analysis of ethylene content and expression of the three ERF genes indicates that ethylene plays a vital role in response to high flooding stress, with ThRAP2.3 being one of the key downstream-response ERF genes regulated by the phytohormone ethylene. Further study is necessary to explore the response mechanism of ThERF15 and ThERF39 to flooding stress and ethylene accumulation, and to verify the roles in flooding stress tolerance of ThERF15, ThERF39 and ThRAP2.3, using genetic transformation, among other strategies.