During growth and development, plants are controlled by many transcriptional factors. The plant-specific homeodomain-leucine (HD-Zip) proteins transcriptional factors affect various biological processes including development, cell division, and responses to abiotic stresses in plants (Elhiti & Stasolla, 2014). HD-Zip proteins are characterized by a C-terminal homeodomain (HD) domain, which is responsible for sequence-specific DNA-binding, and an adjacent leucine zipper (LZ), which facilitates protein–protein interactions. Based on sequence similarity and additional conserved motifs, HD-Zip proteins can be classified into four subfamilies: HD-Zip I, II, III, and IV (Ariel et al., 2007). Compared with HD-Zip I proteins, HD-Zip II proteins contain conserved CPSCE motifs and N-terminal consensus sequences. Both HD-Zip III and HD-Zip IV proteins have a steroidogenic acute regulatory protein (StAR)-related lipid transfer domain (START), which has an unknown function in plants. Compared with HD-Zip III proteins, HD-Zip IV proteins lack the C-terminal MEKHLA motif, but contain the PAS domain, which may be involved in the circadian clock and signal transduction (Mukherjee & Bürglin, 2006; Hefti et al., 2004). To date, many HD-Zip genes have been identified, and many attempts have been made to elucidate the functions of various HD-Zip genes in several species, such as Arabidopsis thaliana (Henriksson, 2005; Nakamura et al., 2006; Ciarbelli et al., 2008), rice (Oryza sativus) (Agalou et al., 2008), maize (Zea mays) (Javelle et al., 2011), tomato (Solanum lycopersicum) (Zhang et al., 2014a; Zhang et al., 2014b), Brassica rapa (Khan et al., 2018), cotton (Gossypium hirsutum) (He et al., 2018), soybean (Glycine max) (Chen et al., 2014), foxtail millet (Setaria italica) (Chai et al., 2018), sunflower (Helianthus annuus) (Manavella et al., 2010; Dezar et al., 2011), poplar (Populus trichocarpa) (Hu et al., 2012), grape (Vitis vinifera) (Jiang et al., 2014), pear (Pyrus betulifolia) (Wang et al., 2015) and peach (Prunus persica) (Zhang et al., 2014a; Zhang et al., 2014b).
Genes within the same subfamily tend to exhibit functional conservation due to structural similarities. Several members of the HD-Zip I subfamily are involved in responses to abiotic stresses and phytohormones, embryogenesis, and de-etiolation and response to blue light. These genes include AtHB6, AtHB7, and AtHB12 from Arabidopsis (Söderman, Mattsson & Engström, 1996; Söderman et al., 1999; Lee et al., 2001), which are induced by drought and abscisic acid (ABA), and HaHB4 from sunflower, which responds to methyl jasmonate (MeJA), ethylene, and drought stress (Manavella et al., 2010). Furthermore, overexpression of AtHB13 in Arabidopsis and HaHB1 in sunflower can significantly improve the cold resistance of transgenic plants (Cabello, Arce & Chan, 2012); the maize HD-Zip gene Zmhdz10 enhances the drought resistance of transgenic rice and Arabidopsis (Zhao et al., 2014). In contrast, MtHB1 reduces the drought resistance of alfalfa by inhibiting lateral root formation (Tina et al., 2017). Moreover, ATHB52 can respond to light and regulate photomorphogenesis and de-etiolation (Henriksson, 2005), and AtHB1 responds to short-day photoperiods and promotes hypocotyl and root elongation (Capella et al., 2015). In addition to the above functions, HD-Zip I genes are associated with the development of plant organs. For example, Arabidopsis hb21/hb40/hb53 triple mutants produced more lateral buds than wild-type plants under short-day conditions, which showed that AtHB21/AtHB40/AtHB53 negatively regulate bud formation (Gonzalez-Grandio et al., 2017).
HD-Zip II genes have key roles in plant responses to light and auxin signalling. For instance, ATHB4 and HAT3 from Arabidopsis, two auxin-induced genes, are involved in shade-induced growth (Bou-Torrent et al., 2012). The sunflower HD-Zip II gene HaHB10 responded to dark and light conditions and was associated with the induction of flowering (Dezar et al., 2011). In addition, there is evidence showing that many HD-Zip II genes are relevant to the regulation of organ growth and development. The AtHAT3, AtHB2, and AtHB4 genes co-regulate shoot apical meristem (SAM) and cotyledon development in Arabidopsis seedlings (Turchi et al., 2013). The Arabidopsis loss-of-function mutant athb4/hat3 results in severely abaxialized leaves (Bou-Torrent et al., 2012).
HD-Zip class III proteins act as regulators of apical meristem development, embryogenesis, lateral organ initiation, and vascular bundle development. The functions of some HD-Zip III genes have been investigated both in herbaceous and woody plants. There are five (REV, PHB, PHV, CAN/ATHB-15, and ATHB-8) and eight (PtrHB1–8) HD-Zip III genes in Arabidopsis and Populus, respectively (Prigge et al., 2005; Hu et al., 2012). REV, PHB and PHV, along with KANADI, are believed to control the abaxial-adaxial patterning of lateral organs, such as leaves and embryos (Emery et al., 2003a). All five members of the HD-Zip III subfamily are thought to be directly involved in vascular development in Arabidopsis by up- or down-regulation. For example, REV, PHB, and PHV promote, while CAN and ATHB-8 inhibit, the formation of interfascicular cambium at the base of the inflorescence stem (Baima, 2001; Prigge et al., 2005). PopREVOLUTA (PRE, PtrHB2), a REV- homologous gene in poplar, plays a key role in regulating the initiation of cambia and the pattern of secondary vascular tissues (Robischon et al., 2011). In contrast, PtrHB5 (PCN), a CAN-homologous gene, inhibits the development of secondary vascular tissue (Du et al., 2011). A PtrHB4 mutation caused a defect in interfascicular cambium in transgenic poplar (Zhu et al., 2018). Because of changes in vascular patterns, some transgenic plants exhibit various morphological characteristics, such as tortuous stems and leaves, dwarfism, and shortened internodes (Du et al., 2011; Zhu et al., 2013; Zhu et al., 2018).
HD-Zip class IV proteins control epidermal cell differentiation and cuticle and trichome formation. Most HD-Zip IV genes are specifically expressed in epidermal or subepidermal layers: these genes include Arabidopsis GL2 (GLABRA2), ATHB-10, ANL2 (ANTHOCYANINLESS 2), ATML1 (ARABIDOPSIS THALIANA MERISTEM LAYER 1), and PDF2 (PROTODERMALFACTOR 2) (Abe, 2003) and OCL1 (OUTER CELL LAYER 1) and OCL4 (OUTER CELL LAYER 4) from maize (Vernoud et al., 2009; Depègefargeix et al., 2011). The atml1/pdf2 double mutant exhibits a few defects in shoot epidermal cell differentiation in Arabidopsis. GL2 and ATHB-10 suppress hair formation to determine trichome and root-hair patterning in A. thaliana by coordinating with single-repeat R3-MYBs (Ohashi et al., 2010). The latest findings suggest that AaHD8 not only promotes the initiation of glandular trichomes but also facilitates leaf cuticle development in Artemisia annua (Yan et al., 2018).
Prunus mume is an economically important fruit and ornamental plant that has been cultivated in China for thousands of years; it is widely distributed and used in East Asia. P. mume has acquired favourable ornamental characteristics, including colourful corollas, a pleasant fragrance, and various flower shapes. In addition, it has different types of plant architecture including the upright, weeping, and tortuous types (Chen, 1996). It is important to improve plant architectural traits that are associated with the development of tissues and organs. HD-Zip transcriptional factors can not only enhance resistance to abiotic stress but also improve morphological features to control the growth and development of multiple tissues and organs, such as the apical meristem, vascular tissue, epidermal tissue, and lateral organs. Compared to abundant annotations of Arabidopsis HD-Zip genes, no P. mume HD-Zip has been analysed in detail. In this study, we identified 32 HD-Zips genes from P. mume and performed detailed structural, evolutionary, miRNA165/166 target gene, and expression pattern analyses in different organs and at different stem development stages. Our study provides an overview of the HD-Zip family in P. mume and information to aid understanding of HD-Zip functions in woody plants.
Materials & Methods
HD-Zip gene identification in P. mume and other species
Candidate genes of the HD-Zip family were extracted from the A. thaliana, P. mume (Zhang et al., 2012), Prunus avium (Shirasawa et al., 2017), and Prunus persica (International Peach Genome et al., 2013) genome databases using the software HMMER (Johnson, Eddy & Portugaly, 2010) based on a Hidden Markov Model (HMM) of the HD-domain profile (PF00046). The genome database of P. mume was downloaded from NCBI (https://www.ncbi.nlm.nih.gov/genome/13911) (Zhang et al., 2012). The genomes of P. avium and P. persica were downloaded from the GDR database (https://www.rosaceae.org/), and the genome of A. thaliana was downloaded from NCBI databases (https://www.ncbi.nlm.nih.gov/). All obtained protein sequences were further analysed by SMART (http://smart.embl-heidelberg.de/) and CD-search (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) softwares to confirm the presence of the leucine-zipper (LZ) domains. Sequences lacking the LZ domains were discarded.
Gene classification and phylogenetic analysis
Multiple sequence alignments of the full-length HD-Zip protein sequences from P. mume, A. thaliana, P. avium, and P. persica were performed by the NAFFT 7 program (Katoh, Rozewicki & Yamada, 2017). The maximum-likelihood (ML) phylogenetic tree was constructed using MEGA 6.0 under the Jones–Taylor–Thornton (JTT) amino acid substitution model with 1,000 bootstrap replicates. The full-length sequences of all HD-Zip proteins used in this study is provided in Data S1.
Gene characteristics, structure, and conserved motif analyses
The genomic sequences, ID numbers, and chromosomal locations of the P. mume HD-Zip genes were obtained from the P. mume genome. The physical parameters of putative proteins, including the lengths of the amino-acid sequences (aa), molecular weights (MWs), theoretical isoelectric points (pIs), and grand averages of hydropathicity (GRAVYs), were calculated using the online ExPASy tool (http://www.expasy.org/tools/protparam.html). The gene structure display server (GSDS) program (http://gsds.cbi.pku.edu.cn/) was applied to illustrate the exon–intron structures by alignment of the CDS with the DNA sequences of individual HD-Zip genes. The orthologues of P. mume HD-Zip in Arabidopsis mapping and GO annotation were conducted using the software Blast2GO (Conesa et al., 2005). The detailed information of all orthologues is provided in Data S2.
The online MEME program was used to display the motif structures of HD-Zip proteins (http://meme-suite.org/tools/meme) using the default parameters, except that the optimum motif widths were set from 6 to 200 residues. SMART (http://smart.embl-heidelberg.de) and CD-search (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) web servers were used for structural motif annotation.
Chromosomal location and gene duplication analyses
The chromosomal location information of HD-Zip genes was obtained from the P. mume genome database, and the chromosomal location map was produced by TBtools (Chen et al., 2018). Tandem duplicated genes were identified, and paralogues were genes separated by five or fewer genes in a 100-kb region, as proposed by Tu et al. (2010). Identification of segmentally duplicated genes was conducted using the Plant Genome Duplication Database (http://chibba.agtec.uga.edu/duplication/index/downloads). Searches for transposons and retrotransposons in DNA sequences 10-kb upstream and downstream of each HD-Zip gene were conducted with RepeatMasker software (http://www.repeatmasker.org/cgi-bin/WEBRepeatMasker/).
Estimation of Ka/Ks ratios
The non-synonymous substitution (Ka) and synonymous substitution (Ks) values of duplicated pairs were calculated on PAL2NAL (http://www.bork.embl.de/pal2nal) (Suyama, Torrents & Bork, 2006; Yang, 2007) after the alignment of protein sequences. The duplication time of duplicated genes was calculated by the following formula: T = Ks∕2λ × 106 Mya (λ = 1.5 × 10−8 for dicots) (Blanc & Wolfe, 2004).
Cis-elements in promoters and miR165/166 target site prediction
The online database PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot, 2002) was employed to investigate putative cis-regulatory elements in the sequences approximately 1,500-bp upstream of the translational start site (ATG).
Full-length HD-Zip nucleic acid sequences (including introns and UTRs) were obtained, and psRNATarget software (http://plantgrn.noble.org/psRNATarget/?function) was used to predict miR165 and miR166 target sites (Wang et al., 2014b).
Expression profile analysis based on transcriptome data
To gain insight into the tissue-specific expression patterns of PmHB genes, our previous RNA-Seq data (GEO No. GSE40162) were used to generate a heat map based on the reads per kilobase per million (RPKM) values of individual gene in all tissue samples. HD-Zip gene expression levels were analysed in five organs of P. mume: roots, stems, leaves, flower buds, and fruits (Zhang et al., 2012). HD-Zip gene expression quantities in flower buds at the endodormancy stage (ED) and natural flush stage (NF) were obtained by RNA sequencing during bud dormancy (Zhang et al., 2018). The heat maps were illustrated using the HemI 1.0 software with the default settings (http://hemi.biocucukoo.org/faq.php).
Preparation of plant materials
P. mume cv. ‘Fei Lve’ was cultivated in the greenhouse of Beijing Forestry University. To investigate the expression of PmHBs, leaf buds and young stems at different developmental stages were sampled from one-year-old branches (Fig. S1). All samples were collected and frozen immediately in liquid nitrogen and then stored at −80 °C until needed for RNA extraction.
RNA extraction and qRT-PCR
Total RNA was extracted from leaf buds and stem tissues using the RNA extraction kit (Tiangen, Beijing, China). After testing for concentration and quality, 500 ng total RNA was used to synthesize first-strand cDNA using the PrimeScript RT Reagent kit (TaKaRa, Dalian, China). The SYBR Premix Ex Taq II kit (TaKaRa, Dalian, China) was used for qRT-PCR according to the manufacturer’s instructions. PP2A and actin were selected as reference genes based on previous reports (Wang et al., 2014a; Xu et al., 2015). Three biological replicates were carried out, and expression levels were calculated by the 2−ΔΔCt method (Schmittgen & Livak, 2008). Gene-specific primers used in the qRT-PCR analyses were designed by IDT (https://sg.idtdna.com/scitools/Applications/RealTimePCR/) and are listed in Table S1.
|Gene name||Gene ID||Locus||Start||End||pI||MW||Strand||Subcellular localization||Gene length (bp)||ORF (bp)||Protein length (aa)||Exons||Introns||Homolog in Arabidopsis|
Identification, characteristics, and phylogenetic analysis in P. mume and other species
In total, 32, 32, 26, and 45 putative proteins containing the HD and LZ domains were identified in P. mume, P. persica, P. avium, and A. thaliana, respectively. The putative proteins of P. mume were named PmHB1–PmHB32 based on their chromosomal positions in the order of successive chromosomes (Table 1). Among the putative PmHBs, physicochemical characteristics were variable. Their encoded proteins had diverse lengths, from 171 aa to 842 aa, with an average length of 463 aa. The molecular weights ranged between 11.98 kDa and 92.32 kDa, and the isoelectric points ranged from 4.64 to 9.31. The number of HD-Zip genes in P. mume was 0.7-fold less than that in A. thaliana (45), whereas it was similar to that in P. persica (32) and less than that in P. avium (26). To investigate the classification of HD-Zip proteins in P. mume and their evolutionary relationships in P. mume, P. persica, P. avium, and A. thaliana, full-length HD-Zip proteins were used to construct a maximum-likelihood phylogenetic tree using MEGA software (Fig. 1A). The HD-Zip proteins of the four species were divided into four subfamilies (HD-Zip I to IV), and the proportions of each subfamily were similar in all four species (Fig. 1B). The HD-Zip I gene family was the largest subfamily in P. mume (37.5%), P. persica (34.4%), and P. avium (38.5%), while HD-Zip IV genes composed the largest proportion in A. thaliana (35.6%). HD-Zip III genes had the fewest representatives in all four species: there were 4 (12.5%), 4 (12.5%), 4 (15.4%), and 5 (11.1%) members in P. mume, P. persica, P. avium, and A. thaliana, respectively.
Chromosomal location and gene duplication of PmHBs
Thirty of the P. mume HD-Zip genes were distributed among all the chromosomes, except for Chr6 (Fig. 2), and only two genes (PmHB31 and PmHB32) remained on unmapped scaffolds. Chr2 and Chr4 contained the most HD-Zip genes, accounting for 18.75% (6) of the entire HD-Zip family; there were five genes on Chr5 and four on Chr7. In addition, Chr1, Chr3, and Chr8 each contained one gene (9.3%).
Because genome duplication is absent in P. mume, segmental duplication, tandem duplication, and transposons are the three mechanisms that resulted in gene expansion in P. mume (Clark & Donoghue, 2018; Del Pozo & Ramirez-Parra, 2015; Qiao et al., 2019). To identify potential segmental duplications in the HD-Zip family, pairs of P. mume and P. persica genes were obtained from the Plant Genome Duplication Database. The distribution of PmHBs relative to their corresponding duplicate genes is illustrated in Fig. 2 and Table 2. Among the 30 PmHBs mapped, seven PmHB gene pairs (43.8%) were found in segmental repeats (PmHB2–PmHB4, PmHB8–PmHB26, PmHB7–PmHB18, PmHB3–PmHB6, PmHB5–PmHB14, PmHB10–PmHB16, and PmHB19–PmHB25), representing all the HD-Zip subfamilies. No genes had been tandemly duplicated, suggesting that tandem duplications were not involved in the expansion of the HD-Zip gene family in P. mume.
|LOCUS 1||LOCUS 2||Ka||Ks||Ka/Ks||Select pressure||MYA|
Interspersed duplication, mainly attributable to transposons, was another mechanism that contributed to the expansion of the HD-Zip gene family. A total of 22 transposons or retrotransposons are found upstream and downstream of 15 PmHB genes. Details about the transposons and retrotransposons are shown in Table S2, implying that transposons may have played an important role in the expansion of HD-Zip genes in P. mume.
The ratio of nonsynonymous to synonymous mutation rates (Ka/Ks) was used to estimate the selection pressure that promoted the evolution of the gene family. The Ka, Ks, and Ka/Ks of 19 orthologous gene pairs between P. mume and P. persica and seven paralogous gene pairs in P. mume were calculated (Table 2). All the Ka/Ks of the paralogous gene pairs were <1, which indicates that the HD-Zip genes of P. mume have mainly experienced purifying selection pressure with limited functional divergence. In addition, most of the Ka/Ks ratios of the orthologous genes were <1, while the Ka/Ks value of PmHB9/ppa016510 m was >1, suggesting that HD-Zip genes evolved primarily under purifying selection that occurred after the divergence of P. mume and P. persica; adaptive evolution may have occurred in PmHB9.
Conserved domains and gene structure analysis of PmHBs
The conserved motifs of the HD-Zip protein sequences were analysed by MEME (Fig. 3A). Motif 1 and motif 2 were found to encode the homeobox domain. Motif 1 was found in all 32 genes. In addition, motif 7 was found to encode the LZ domain. Motifs 3, 4, 9, and 12 encode START domains and are present in subfamilies III and IV. Motif 8 was only found in subfamily III and encodes the MEKHLA domain. In addition, it was suggested that some unannotated subfamily-specific motifs might contribute to the diversification of the subfamily. A schematic representation of the HD-Zip family based on the HD domain, LZ domain, START domain and MEKHLA domain regions is shown in Fig. 3B.
The number of introns is similar in the different subfamilies (Fig. 3C). Most of the HD-Zip I genes have two introns, except PmHB12. In the HD-Zip II subfamily, PmHB3, 6, 9, and 31 have two introns, while PmHB19, 25, 28, and 32 have three introns. HD-Zip IV subfamily genes have 8–11 introns, whereas all the HD-Zip III genes have 17 introns, which was the maximum number among the HD-Zip family in P. mume. The results above divided the HD-Zip genes into four subfamilies, which were consistent with the result of the phylogenetic analysis.
MicroRNA165/166 target-site prediction and cis-acting elements of the PmHBs promoters
Previous studies revealed that the activity of HD-Zip III proteins was regulated by miR165/166 (Kidner & Martienssen, 2004; Mallory et al., 2014). Here, miR165 and miR166a–e target sites were found in four P. mume HD-Zip III subfamily genes (PmHB1, PmHB5, PmHB11, and PmHB14) (Fig. 4), suggesting that P. mume HD-Zip III subfamily genes may be regulated by miR165 and miR166.
The cis-acting elements of promoters are essential for determining tissue-specific expression and are involved in the regulation of gene expression. A region of 1500 bp upstream of the start codon (ATG) of HD-Zip genes was used to analyse cis-acting regulatory elements. Many of the cis-acting elements identified in the PmHB gene promoters were related to the light response, including the AE-box, Box 4, the I-box, the GATA-motif, the G-box, the GT1-motif, and the TCT-motif (Fig. 5), suggesting that PmHBs may respond to light and participate in plant light morphogenesis. Because many cis-elements associated with hormone responses (the P-box, TATC-box, TCA-element, ABRE, TGA-element, AuxRR-core, TGACG-motif, CGTCA-motif, and ERF) were distributed in the promoter regions, the expression of PmHB genes may be induced by several plant hormones, such as salicylic acid, MeJA, abscisic acid, gibberellins, and auxin. In addition, the presence of shoot- and meristem-specific expression elements (as-2-box, CAT-box, and O2-site) in the promoters of several PmHB genes indicated that those genes might participate in regulating the development of the meristem and stem. HD-Zips are deemed to be crucial TFs involved in meristem and stem development.
Expression analysis of PmHBs in different tissues of P. mume
To explore the expression patterns of PmHBs in different organs and their biological functions during plant growth and development, the expression levels of PmHBs in five tissues (leaf, flower bud, fruit, root, and stem) were analysed by RNA-seq. Heatmaps showed that HD-Zip gene expression levels are diverse in the five tissues (Fig. 6). PmHB29 transcripts were not present in fruits and roots, and PmHB9 was not present in buds and stems; however, other HD-Zip genes could be detected in all five tissues. These results indicated that some HD-Zip genes have multiple functions in the development of different organs. PmHB2 and PmHB26 were highly expressed in all tissues, whereas PmHB9, 13, and 27 were expressed at low levels in all tissues. Twelve genes (PmHB1, 5, 7, 11, 14, 16, 18, 20, 22, 25, 26, and 27) were most highly expressed in the stem, which showed that these genes might have functions in regulating stem development. It is notable that PmHB10, 16, 20, and 24 were expressed in all aboveground organs but were expressed either at low levels or not at all in roots. In addition, some genes of the same subfamilies exhibited similar expression patterns. For instance, all HD-Zip III subfamily genes (PmHB1, 5, 11, and 14) showed high RPKM values in the stem, and five HD-Zip IV subfamily genes (PmHB16, 17, 20, 22, and 24) showed high RPKM values in both leaf and stem tissues.
HD-Zip genes may also regulate flower bud dormancy. Screening for differentially expressed genes (DEGs) by RNA sequencing in flower buds at the endodormancy (ED) and natural flush (NF) stages was carried out in our previous study (Zhang et al., 2018). Eleven genes (PmHB6, 7, 8, 11, 19, 26, 27, 29, 30, and 31) were more highly expressed at the ED stage, and four genes (PmHB3, 17, 18, and 23) were upregulated at the NF stage (Table S3).
Expression pattern analysis of PmHBs in stem development
To verify the functions of PmHBs in the development of leaf buds and stems, RNA-seq on leaf buds and stem tips was conducted in P. mume. HD-Zip family genes that were differentially expressed the between leaf bud and stem tip are shown in Table S4. Seven genes (PmHB10, 16, 20, 22, 24, 27, and 29) were upregulated in leaf buds, while another six genes (PmHB3, 13, 17, 18, 21, and 30) were upregulated in stem tips. The results above suggest that HD-Zip genes play putative roles in the development of leaf buds and stems.
Expression levels of 10 selected genes in different organs (leaf, leaf bud, and differentiating stem) were determined by qRT-PCR (Fig. 7). The results showed that genes in the HD-Zip III subfamily had similar expression patterns, which were highly expressed in both LBSII and stems, and these patterns were not significantly different at different stages of stem development; PmHB7 and PmHB18, two HD-Zip I genes, were predominantly expressed in S6–10 in stem, while another HD-Zip I gene, PmHB17, exhibited high levels of expression in leaf and in S11–15. PmHB20 of the HD-Zip IV subfamily had higher expression levels in leaves, leaf buds, and large leaf buds and lower expression levels in different stages of stem development. PmHB19 is a HD-Zip II gene with high expression levels in both large leaf buds and stem tips (S1–5). Another HD-Zip II gene, PmHB25, was expressed maximally in S11–15, while it exhibited its lowest expression levels in LBSI and LBSII (Fig. 7).
Expression in different tissues of the stem (xylem, cambium, and bark) were determined by qRT-PCR (Fig. 8). HD-Zip III subfamily genes (PmHB1, 5, 11, and 14) may participate in xylem differentiation because they are expressed predominantly in the xylem. However, PmHB7, 17, 18, 19, 20, and 25 were mainly expressed in bark, suggesting that they might be involved in the development of phloem, endodermis, or epidermis (Fig. 8).
Interaction network between ten P. mume and A. thaliana HD-Zip proteins
The GO annotation of HD-Zip genes in P. mume were conducted by Blast2GO software basing on the orthologues in Arabidopsis (Fig. S2). PmHBs were divided into three GO categories: biological processes, cell components and molecular functions. In the biological process category, biological regulation was the most highly represented groups. The cell, cell part, and organelle classes were abundant in cellular component category. Within the molecular function category, genes that corresponded to binding were the most abundant.
To further predict P. mume HD-Zip proteins functions during plant stem development, the protein–protein interactions between ten HD-Zip proteins of P. mume and A. thaliana were analysed using STRING software (http://string-db.org) (Fig. 9). The homologous genes were matched with the highest bit score, the results were corresponding to those of Blast2GO and the lines of different colours represent different types of evidence for the proteins interaction. Ten HD-Zip proteins were found to interact in a network, suggesting that they might have complex connections in the co-regulation of stem development. Previous studies showed that HD-Zip III activity can promote the expression of ZPR proteins in A. thaliana (Husbands et al., 2016; Wenkel et al., 2007). Specifically, REV and PHB can interact with ZPR3 to participate in the regulation of leaf polarity and vascular development (Kim et al., 2008; Wenkel et al., 2007). Moreover, interactions between the HD-Zip III and KANDAI proteins to regulate vascular development have been reported in A. thaliana (Emery et al., 2003a; Emery et al., 2003b). The HD-Zip III subfamily genes PmHB1, 5, 11, and 14 had highly similar amino sequences to those of REV, HB-8, PHB, and ATHB-15, respectively, and might regulate leaf and vascular development. In addition, HB-2 was induced by the PIFs and participated in the shade avoidance response (Kunihiro et al., 2011). PmHB25, which is closely related to HB-2, may play an important role in regulating light-mediated stem elongation.
HD-Zip, containing a homeodomain (HD) domain and an adjacent leucine zipper (LZ), is a major group of transcriptional factors present in higher plants that has a significant impact on the light reduction, abiotic stress, and organic development processes of plants. Most studies on HD-Zip proteins have focused on herbaceous plants, while the functions of only a few HD-Zip proteins have been verified in woody plants, such as poplar and peach (Hu et al., 2012; Zhang et al., 2014a; Zhang et al., 2014b). However, a great deal is unknown about these genes in P. mume. A total of 32 genes were identified, which is 0.7-fold more than the number found in Arabidopsis (Emery et al., 2003b; Henriksson, 2005; Nakamura et al., 2006; Ciarbelli et al., 2008), 0.5-fold more compared to Populus (Hu et al., 2012), and similar to the number found in P. persica (Zhang et al., 2014a; Zhang et al., 2014b). The number of HD-Zip genes varies widely among these species. Nevertheless, the P. mume genome (280 Mb) is approximately two-fold larger than that of Arabidopsis (125 Mb) (Zhang et al., 2012; Initiative, 2000). To explore the phylogenetic relationship between the HD-Zip proteins of P. mume and Arabidopsis, PmHB orthologues in Arabidopsis with the highest levels of homology were extracted (Table 1), showing that a single HD-Zip gene in Arabidopsis had multiple orthologues in P. mume. For instance, ATHB-6, HDG11, ANL2, and HOX11 had two orthologues each in P. mume, whereas HAT22 and HAT3 had three orthologues each in P. mume. These results suggested that, such as peach, some HD-Zip genes may have been lost in the P. mume lineage during evolution (Zhang et al., 2014a; Zhang et al., 2014b). The number of HD-Zip genes in P. mume was similar to that in peach, which might be because of their closer genetic relationship.
The 32 putative genes were classified into four subfamilies according to the results of phylogenetic analysis, and the length and intron number of the genes in different subfamilies were conserved (Figs. 1 and 3). For example, the lengths of HD-Zip I and II genes were shorter than those of HD-Zip III and IV. HD-Zip III genes have 17 introns. In addition, proteins in the same subfamily showed similar structures. In general, most of the motifs and domains in the same subfamily were similar, which may be consistent with the phylogenetic analysis. However, some motifs were specific to subfamilies. For example, motif 18 was only detected in the HD-Zip II subfamily; motifs 3 and 4 were found in the HD-Zip III and IV subfamilies, and motifs 5, 6, 9, 10, 12, 15, and 17 were only found in the HD-Zip IV subfamily, which might be associated with the functional diversity of the HD-Zip family. These findings clearly support the results of phylogenetic analysis. Tandem duplication, segmental duplication, and transposons have resulted in the genome evolution. It was revealed that segmental duplication and transposons may have contributed to the expansion of the HD-Zip gene family in P. mume, accounting for 92% of duplications (Fig. 2 and Table 2).
Proteins can interact with other proteins to co-regulate biological processes. A network of ten P. mume and Arabidopsis HD-Zip TF pairs was established. In response to shade, plants can react by altering elongation during growth and development. It has been found that HD-Zip genes positively regulate the plant response to shade. In Arabidopsis, the HD-Zip II gene HB-2 can be activated by PIF4 and PIF5, two phytochrome-interacting transcriptional factors that participate in the regulation of hypocotyl elongation under short-day conditions (Kunihiro et al., 2011). As the orthologues of HB-2 and HAT3, PmHB25 and PmHB19 may also play roles in shade avoidance. Moreover, Arabidopsis hat3/athb4/athb2 loss-of-function triple mutants showed developmental defects in white light, such as adaxialized leaves and reduced shoot apical meristem activity (Turchi et al., 2013). Previous research showed that the expression of several HD-Zip II transcriptional factors, including HAT3 and HB-2/HAT4, was regulated by the HD-Zip III gene REVOLUTA (REV) (Bou-Torrent et al., 2012). Mutation of REV, PHABULOSA (PHB), and PHAVOLUTA (PHV) enhanced the bilateral symmetry and SAM defects in hat3/athb4 mutants. These results indicated that HD-Zip II and HD-Zip III genes cooperate in establishing bilateral symmetry and patterning in the embryo, as well as in controlling SAM activity in plants. Six genes in both the HD-Zip II and III subfamilies (PmHB1, 5, 11, 14, 19, and 25) had higher expression levels in leaf buds before leaf unfolding (LBSII) or in stem tips (S1–5), suggesting that PmHB1, 5, 11, 14, 19, and 25 might participate in the formation and maintenance of SAM polarity in P. mume.
HD-Zip III transcriptional factors also play roles in controlling primary and secondary vascular cell differentiation and SAM (shoot apical meristem) maintenance (Baima, 2001; Kim et al., 2005), which is essential for the formation of plant architecture. Overexpression of ATHB-8 promotes vascular cell differentiation and xylem tissue production in inflorescence stems of Arabidopsis (Baima, 2001). HD-Zip III TFs are key factors that regulate the secondary growth of vascular bundles in Populus. PtrHB7, an ATHB-8 homologue, regulates the balance of vascular patterning by inhibiting the differentiation of secondary xylem while promoting the differentiation of secondary phloem (Zhu et al., 2013). REV, which has overlapping functions with PHV and PHB, regulates meristem initiation in lateral organs in opposition to KANADI (Ilegems et al., 2010). Overexpression of PRE, the Populus homologue of the Arabidopsis REV gene, in Populus causes the patterning of secondary vascular tissues to be altered (Robischon et al., 2011). In addition, genes of the four HD-Zip III subfamilies of P. mume were expressed preferentially in xylem and cambium; thus, we speculated that four genes from the HD-Zip III subfamily were closely associated with stem development in P. mume. The expression of HD-Zip III genes was controlled by multiple molecular mechanisms. First, the small ZIP protein ZPR3 reduced HD-Zip III protein activity by direct interaction to form nonfunctional heterodimers in the SAM of Arabidopsis (Kim et al., 2008). Moreover, overexpression of KANADI negatively affected the expression of HD-Zip III genes during embryogenesis, whereas HD-Zip III genes and KANADI co-regulated auxin flow by controlling the distribution of PIN-FORMED1 (PIN1) (Izhaki & Bowman, 2007). In addition, protein-interaction networks showed that REV, HB-8, PHB, and ATHB-15, the orthologues of PmHB1, PmHB5, PmHB11, and PmHB14, respectively, interacted with ZPR3, and ATHB-15 interacted with KANADI (Fig. 9). Finally, MiR165/166 negatively regulated the activity of HD-Zip III genes in both Arabidopsis (Kim et al., 2005; Mallory et al., 2014) and Populus (Ko, Prassinos & Han, 2006). Because the target sites were present in the sequences of HD-Zip III genes, we speculated that HD-Zip III genes might also be regulated by MiR165/166. However, the molecular mechanisms of the HD-Zip III gene regulation of vascular development and interactions between proteins need further investigation.
HD-Zip IV proteins have been reported to regulate the differentiation and maintenance of outer cell layers (Nakamura et al., 2006; Yan et al., 2018). PDF2 and ATML1 regulate shoot epidermal cell differentiation and embryo development in Arabidopsis (Ogawa et al., 2015). In P. mume, the expression level of PmHB20 in bark was more than 2,000 times that in cambium and xylem (Fig. 8), which is similar to Arabidopsis HD-Zip IV genes that are preferentially or specifically expressed in plant epidermal or sub-epidermal cells (Rombolácaldentey et al., 2014).
Thirty-two HD-Zip genes were identified in P. mume. Phylogenetic analysis showed that there were four subfamilies in the HD-Zip family: subfamilies I, II, III, and IV. Gene structure and motif analysis suggested that genes in the same subfamily were conserved. Furthermore, expression profile analysis showed that HD-Zip genes may be involved in stem development and flower bud dormancy. Otherwise, all HD-Zip III genes were highly expressed in cambium and xylem, suggesting that HD-Zip III genes may be indispensable in vascular differentiation. The prediction of microRNA target sites indicated that HD-Zip III subfamily genes may be regulated by miR165/166. Promoter analysis showed that HD-Zip III genes may be involved in responses to light, hormones, and abiotic stressors, as well as plant development. Our study identified the members of the HD-Zip family in P. mume and laid the foundation for molecular breeding programmes for woody ornamental plants.
Sequence of forward and reverse primer pairs for real-time quantification PCR
Transposons and retrotransposons in genomic sequences of 10-kb upstream and downstream of each HD-ZIP gene
Differentially expressed genes (DEGs) in flower buds with endodormancy stage (ED) and natural flush stage (NF) in P. mume
Differentially expressed genes (DEGs) in leaf bud and stem tip of P. mume
Samples of leaf buds and young stems in different stages
L, leaf; LBSI (leaf buds in stage I), leaf bud in sprouted phase; LBSII (Leaf bud in stage II), leaf buds before leaf unfolding; S1-5, stems from 1 to 5 internodes; S6-10, stems from 6 to 10 internodes; S11-15, stems from 11 to 15 internodes