Genome-wide identification, characterization and expression analysis of the HD-Zip gene family in the stem development of the woody plant Prunus mume

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λ × 10⁶ Mya (λ = 1.5 × 10⁻⁸ 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.

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.

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.