To identify members of the TIFY gene family in poplar, the Hidden Markov Model (HMM) profile of the TIFY domain (PF06200) was downloaded from the pfam database (Pfam²). The reference genome P. trichocarpa was downloaded from the Phytozome database³. Genome annotation version 3.0 of P. trichocarpa was used in this study. To search the 73,013 predicted proteins of the P. trichocarpa genome annotation version 3.0, 759 seed sequences in the PF06200 HMM profile were used. Two softwares, blastp and HMMER, were used to perform a search for proteins harboring the TIFY domain. Proteins of P. trichocarpa that showed E-values above 1e-6 in the search results of blastp or HMMER were considered to be candidate TIFY domain genes. Of the proteins from the same gene model, only the longest ones (the “0.1” gene model of TIFY genes in the P. trichocarpa genome annotation version 3.0) were kept for further study. To confirm the existence of the TIFY domain and the Jas, CCT and ZML motifs, candidate TIFY domain proteins were used to search the Pfam database. The Pfam accession numbers for the three motifs are PF09425, PF06203 and PF00320, respectively. We also checked for the existence of the PPD domain/motif in each of the candidate TIFY genes. We noted that there is no sequence information for the PPD domain/motif in the Pfam database. The previously annotated PPD genes in Arabidopsis and grape were used to build the conserved PPD sequence. Briefly, the previously annotated PPD sequences were first aligned by ClustX2.1, and an msf result was produced and imported into hmmbuild implemented in HMMER. The resulting PPD HMM profile was used to search all of the candidate TIFY genes.

The software MEGA5 was used to analyze the phylogeny of all of the TIFY genes (Tamura et al., 2011). The protein sequences of all of the identified poplar TIFY genes were imported into ClustalX2.1 to perform a complete alignment, and the resultant multiple alignment file was imported into MEGA5. Unrooted phylogenetic trees were constructed using the Neighbor-Joining (NJ) method, and the bootstrap test was carried out with 1,000 iterations. For constructing phylogenetic trees among the four species, the protein sequences of the three other species were obtained from Phytozome database⁴. The deduced TIFY protein information was obtained from Arabidopsis thaliana (Bai et al., 2011), Oryza sativa (Bai et al., 2011) and Vitis vinifera (Zhang et al., 2012).

Gene structures and their chromosomal locations were obtained from the P. trichocarpa genome annotation version 3.0. The gene structures were displayed using a custom R script. The duplication pattern for all of the poplar genes was analyzed using MCScanX software. Briefly, the 41,335 poplar gene models (the “0.1” gene model in the P. trichocarpa genome annotation version 3.0 represents all genes models if a gene has more than one alternative transcript) were extracted from the poplar genome. An all-vs.-all local blast for the 41,335 poplar gene models using the local blast software with E-values under 1e-4 was carried out. The blast output was imported into MCScanX software following the procedure. All of the 41,335 genes were classified into four types using default criteria, including segmental, tandem, proximal, and dispersed duplications. The duplication pattern for each poplar TIFY gene was then determined according to the results.

The PopGenExpress database was employed to investigate the expression levels of all of the poplar TIFY genes in different tissues⁵ (Wilkins et al., 2009). Briefly, the correspondence between microarray probes and poplar genes was obtained from the annotation files deposited for the microarray data. Because the annotation released for the poplar microarray used in the PopGenExpress database was based on P. trichocarpa genome annotation version 2.0, the relationship between the probes and poplar genes annotated in genome version 3.0 was transformed based on the correspondence between versions 2.0 and 3.0 of the P. trichocarpa genome annotation. For genes that had no corresponding probes in genome version 3.0, an additional reciprocal blast with released sequences for all of the probes was employed. If a given gene and probe pair were each other’s top hits in the reciprocal blast results, then the probe was considered to represent the gene. The relative expression levels of all of the TIFY genes in xylem, roots, mature and young leaves, male and female catkins were obtained from the PopGenExpress database. A heat map for the expression levels of all of the TIFY genes in the selected tissues was prepared with Cluster 3.0 software (de Hoon et al., 2004); the dendrogram was visualized using JavaTreeview⁶.

The hybrid poplar variety “NL895” (P. euramericana) was used in this study. The female parent of NL895 is P. deltoids cv. Lux” (I-69) and the male parent is also a hybrid variety “I-45” (P. × euramericana). Tissue culture seedlings of “NL895” were grown in 3-L pots containing a sand-peat (50:50, v/v) mixture, and the pots were placed in a greenhouse under a 16/8 h photoperiod. The light intensity in the greenhouse was set at 12,000 Lx. The temperatures were set at 28 and 25°C for day and night, respectively. The seedlings were watered twice per month with Hoagland’s complete nutrient solution. Additionally, the seedlings were watered with distilled water 2–3 times a week. Seedlings with 6–10 fully expanded leaves were used for the treatments.

For the plant hormone treatments, “NL895” seedlings were treated with 0.2 mM jasmonic acid (JA) and 0.5 mM salicylic acid (SA). Leaves with LPIs (leaf plastochron indices) ranging from 4 to 10 were sprayed with JA or SA solution. To ensure that the hormone treatments were uniform, leaves were sprayed on both sides just until the leaf surfaces began to form small drops. Leaf samples were collected at 0, 2, 6, 12, and 24 h after treatment. Sprayed leaves from each individual poplar plant were pooled into one biological sample. Three replicates were carried out for each treatment.

For the leaf rust pathogen infection treatments, a virulent M. larici-populina strain used as fungal material in this study was obtained from P. simonii Carr., which showed typical leaf rust symptoms (yellow bubbles on the lower surface of leaf) in the summer of 2015. This strain was also used as the fungal material in another study (Wang et al., 2017). The urediniospores of M. larici-populina isolates were allowed to multiply on 1-year-old potted P. simonii Carr. cuttings and diluted into an urediniospore suspension in agar-water (0.1% agar in distilled water). The concentrations of the urediniospore suspensions were 1–2 mg/ml. Fully expanded leaves from “NL895” with leaf plastochron indices (LPIs) of 4–10 from each branch were spray-inoculated on their abaxial surfaces with the prepared urediniospore suspension. Leaves from one separate cutting were considered to be one sample. Leaf samples were collected after 0, 2, 4, and 8 days treatments. The samples were immediately frozen in liquid nitrogen and stored at -80°C until required for further analysis.

Total RNA was isolated using an RNAprep Pure Kit (for Plants) according to the manufacturer’s protocol [TIANGEN Biotech (Beijing) Co., Ltd., Beijing, China]. The quality of all of the RNA samples was examined by performing agarose gel electrophoresis. The first strand cDNA was synthesized using the TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix Kit (TransBionovo Co., Ltd., Beijing, China). The primers used for the qRT-PCR assay of TIFY genes are listed in Supplementary Table S2. PCR products from cDNA or genomic DNA templates using these primer pairs were initially sent for sequencing to confirm that the correct targets had been amplified (data not shown). The correct primer pairs were used for subsequent qRT-PCR assays that were performed on the LightCycler 96 (Roche) platform by using the FastStart Essential DNA Green Master Mix (Roche). For each sample, two reference genes (ACTIN and 18S) were used to standardize the mRNA abundance, and three replications were performed. The 2^-ΔΔCq method was used to calculate the relative gene expression based on the qRT-PCR data (Livak and Schmittgen, 2001).

A genome-wide search of TIFY family genes in poplar yielded a total of 24 non-redundant genes that were found to harbor the TIFY domain. The 24 non-redundant genes represent a total of 99 transcripts in the P. trichocarpa genome annotation version 3.0 (Table 1) and all the 99 transcripts have the TIFY domain. To avoid complexity in the subsequent study, only the “0.1” gene model for the 24 non-redundant genes were kept in further analyses. These genes were considered to be TIFY family genes (Table 1). Their phylogeny and gene structures are shown in Figures 1A,B, respectively. Of these 24 TIFY domain-containing genes, 12 of the predicted proteins contain both the TIFY domain and the Jas motif; these proteins were designated as PtJAZ1 through 12. Seven proteins have both a TIFY domain and CCT and ZML motifs, and one protein has both a TIFY domain and a CCT motif; thus, these eight proteins were designated as PtZML1 through 8. Two proteins contain the TIFY domain and the PPD motif, and they were designated as PtPPD1 and PtPPD2. Note that PtPPD1 also has a Jas motif in its protein structure. The remaining two proteins contain only the TIFY domain, and they were designated as PtTIFY1 and PtTIFY2. All of the domains or motifs located in the 24 TIFY genes are illustrated in Figure 1C. The nomenclatures of these genes were designated based on both their subfamilies and their locations on the poplar chromosomes (Figure 2). The correspondence of our nomenclatures for the identified TIFY family genes and their original gene IDs released by the P. trichocarpa genome annotation version 3.0 are shown in Table 1. In previous report, 25 TIFY proteins were identified using the P. trichocarpa genome version 1.0 annotations when global comparisons of the gene family were performed in a number of plant species (Bai et al., 2011). When comparing our results with theirs, we found PtZML4 and 5 were corresponding to 3 ZML proteins in Bai‘s results; while the other 22 TIFY genes are consistent in the two identifications.

According to Figure 1B, all of the TIFY family genes have more than one exon. PtJAZ9 has the smallest number of exons at 2, whereas PtZML5 has the largest number of exons at 11. The genomic sequences of most of these genes are less than 6 kb in length, but PtZML4 has the longest genomic sequence of approximately 12 kb. The gene PtJAZ9 has the shortest open reading frame (ORF) of 149 AA, whereas PtPPD2 has the longest ORF of 454 AA. Insights into the phylogeny and gene structures of genes has indicated that proteins that are closely related phylogenetically tend to have similar gene structures. For example, PtJAZ1 and PtJAZ4 have the same number of exons in the same arrangement and tend to cluster together in the phylogenetic tree; this phenomenon was also observed for PtTIFY1 and PtTIFY2 and PtZML1 and PtZML3. In order to further insight into the phylogeny of TIFY proteins within plant, TIFY family proteins of three other representative species, including Arabidopsis, rice and grapes, were also employed to construct a phytogenic tree. Clearly, TIFY proteins in the four plants can be grouped into four clades (Supplementary Figure S1). According to this figure, proteins in the clade 1 and 2 tend to only have TIFY domain and have ZML motif, respectively; while proteins in the clade 3 and 4 have no clear difference when comparing with their domains or motifs.

The chromosomal locations of these 24 TIFY family genes is shown in Figure 2 and Table 1, and they are found on 13 of the 19 chromosomes. Chromosomes 4, 9, 13, 14, 16, and 19 contain no TIFY genes. Because duplication usually contributes to the expansion of gene families, we investigated the duplication patterns of each of the genes. In total, 19 genes were produced by whole genome duplication (WGD), and two genes were produced by tandem duplication. The other three genes did not show duplication patterns. The two tandemly duplicated genes, PtZML2 and Pt ZML6, were found to be tandemly copied from PtZML1 and PtZML5 (Figure 2), respectively, whereas both PtZML1 and PtZML5 were produced by WGD. This result suggested that PtZML1 and PtZML5 were the ancestors of PtZML2 and PtZML6, respectively. The gene duplication patterns for all 24 of the TIFY genes showed that WGD played large roles in this gene family in poplar.

Since the transcriptional abundance of a gene in different tissues is usually indicative of its function, we investigated the expression profiles of the 24 TIFY family genes in poplar. By analyzing of the publicly available expression data from the PopGenExpress database, we were able to obtain the relative expression levels of all 24 TIFY family genes in six tissues: xylem, roots, mature and young leaves and male and female catkins (Supplementary Table S1 and Figure 3). All 24 genes showed very high gene expression levels in female (7.11 ± 4.54) and male catkins (5.84 ± 1.63), whereas the lowest expression levels were found in roots (Supplementary Table S1). This result may suggest that the TIFY family of genes plays important roles in poplar catkin development.

Genes that are expressed abundantly in one or several specific tissues usually indicate their function related to the formation of these organs or the corresponding plant development. Further insights into the expression levels for each individual gene indicated that TIFY genes showed extremely high expression levels in some tissues. We filtered these types of genes when their expression was greater than two-fold the average expression level of all TIFY genes in the tissue. In total, seven TIFY genes were preferentially expressed in a given tissue, and some tissues had more than one preferentially expressed gene. PtJAZ6 and PtJAZ9 were preferentially expressed in xylem; PtJAZ2 and PtJAZ9 were preferentially expressed in roots and mature leaves; PtJAZ1 and PtJAZ2 were preferentially expressed in both male and female catkins, and PtJAZ9 was preferentially expressed in female catkins; PtJAZ7 and PtJAZ11 were preferentially expressed in young leaves (Supplementary Table S1 and Figure 3). These data suggest that PtJAZ1, PtJAZ2, PtJAZ6, PtJAZ7, and PtJAZ9 could be responsible for the formation of these tissues or contribute to the functions of these tissues. By summarizing these results, we were able to find that a large proportion of genes in the JAZ subfamily were filtered out by our method, although some genes in other subfamilies also showed relatively high abundance in the seven selected tissues. This suggests that most of the JAZ subfamily in the TIFY gene family plays important roles in poplar development. In addition, the expression values for the two reference genes, Actin and 18S, across the five developmental tissues were also obtained and listed in Supplementary Table S1. The data revealed these two genes were expressed similarly in different tissues and it suggested that these two genes were suitable for the subsequent qRT-PCR assay.

JA and SA are hormones that play important roles in signal transduction when plants are challenged with biotic and abiotic stresses. Treating plants with JA, SA or their derivatives MeJA and MeSA can modulate the symptoms of pathogenic or herbivory damage on plants (Yamada et al., 2012; Zhang et al., 2016). To investigate which genes in TIFY gene family can respond to these two plant hormones, we examined the expression profiles of TIFY genes in poplar leaves after treatment with MeJA and SA solutions. The relative expression levels were assayed by qRT-PCR. In the description to convenience, CK indicated untreated samples; JA-2 h, JA-6 h, JA-12 h, JA-24 h indicated samples collected from leaves after 2, 6, 12 and 24 h of MeJA treatment, respectively; and SA-2 h, SA-6 h, SA-12 h, and SA-24 h indicated samples collected from leaves after 2, 6, 12, and 24 h of SA treatment, respectively. In total, we were able to examine the expression profiles for 20 of the 24 TIFY genes; the expression profiles of PtJAZ8, PtPPD2, PtZML4 and PtZML5 were not examined due to unsuccessful primer design. The relative expression profiles of all 20 genes for all nine samples are shown in Figure 4.

Fold-changes in expression larger than 2.0 or less than 0.5 and also with a p-value less than 0.05 between CK and time points after hormone treatment was considered that gene expression was influenced by MeJA or SA treatment. According to the data presented in Figure 4 and this criterion, gene expression patterns could be classified in four groups. The first group showed that both MeJA and SA increased gene expression at some time points but decreased it at other time points. This group included PtJAZ1, PtJAZ4 and PtJAZ7. The second group showed that both MeJA and SA decreased gene expression at all of the time points (except SA-24 h). This group included PtJAZ2, PtJAZ9, PtZML3, PtZML6, PtZML8, PtTIFY1 and PtTIFY2. Expression of most of the genes in this group was significantly decreased at JA-2 h, JA-6 h, JA-12 h, JA-24 h, SA-2 h, SA-6 h, and SA-12 h. In the third group, MeJA decreased gene expression at most time points, whereas SA increased it at some time points and decreased it at others. This group included PtJAZ3, PtJAZ5, PtJAZ10, PtPPD1, PtZML1, PtZML2 and PtZML7. In the fourth group, MeJA increased gene expression at some time points and decreased gene expression at the other time points, whereas SA decreased gene expression at all of the time points. This group included PtJAZ6, PtJAZ11 and PtJAZ12.

Previous reports suggested JAZ protein might act as TFs in response to rust fungi (Duplessis et al., 2009). To identify TIFY genes that respond to poplar leaf rust disease, we collected leaf samples from the poplar hybrid variety “NL895” that had been infected with M. larici-populina (see Materials and Methods); these samples were designated as CK, PI2, PI4, and PI8. The CK represented samples collected immediately after leaf rust infection, while PI2, PI4, and PI8 represented samples collected from leaves treated by urediniospore suspension of M. larici-populina after 2, 4, and 8 days, respectively. The symptom of leaf rust on leaves of “NL895” at these four selected time points were different. No visible symptom was observed at CK and 2 dpi; while very small white milk-white and plenty of uredinia could be observed at 4 and 8 dpi, respectively (data not shown). Previously, 48 hpi (equal to 2 dpi) was reported to be a key time point for the rust development (Rinaldi et al., 2007). Therefore, our plant material could be able to represent all key developmental stages of M. larici-populina on poplar leaves.

The relative expression of genes for PI2, PI4, and PI8 versus CK for all of the 24 TIFY genes was determined by qRT-PCR assays. We were able to test the expression of 18 of the 24 TIFY genes (Figure 5). According to the expression profiles of these 18 TIFY genes and using the criterion of fold changes less than 0.5 or larger than 2 and also with a p-value less than 0.05 when comparing CK with the other three treatments, we were able to classify them into four groups. The first group showed that gene expression was reduced by M. larici-populina infection at all the three stages of PI2, PI4, and PI8 when comparing with CK. This group included PtJAZ5, PtJAZ7, PtZML3, PtZML5 and PtZML6. The second group showed similar gene expression level of CK and PI8, while gene expression was largely reduced at PI2 and PI4 stages. This group had five genes and it included PtJAZ1, PtJAZ11, PtPPD1, PtTIFY1 and PtZML2. The third group showed that gene expression was similar at CK, PI2 and PI4 and highly increased at PI8. This group had six genes and it included PtJAZ2, PtJAZ6, PtJAZ9, PtJAZ12, PtTIFY2 and PtZML1. The fourth group showed that gene expression was reduced at PI2 and PI4 and highly increased at PI8. This group had two genes and it included PtJAZ3 and PtZML8. Based on these data, we could conclude that all tested TIFY genes respond to M. larici-populina infection and they also showed different responsive patterns.