Genome-wide identification and characterization of TIFY family genes in Moso Bamboo (Phyllostachys edulis) and expression profiling analysis under dehydration and cold stresses

Data resources

The whole genome dataset, full length cDNA, and expressed sequence tags (ESTs) of P. edulis were downloaded from the bamboo genome database (www.bamboogdb.org/) or National Center for Biotechnology Information (www.ncbi.nlm.nih.gov/). These include 31,987 protein-coding genes predicted from whole genome sequences, 10,608 low redundant full-length cDNA sequences and 38,000 ESTs from leaf, shoot, and seedling libraries. The deduced amino acid sequences of TIFY genes in A. thaliana genome and rice genome (O. sativa subsp. indica) were according to Vanholme et al. (2007) and Ye et al. (2009), respectively. The corresponding nucleotide sequences of coding region (cds) and genomic sequences were obtained from The Arabidopsis Information Resource (TAIR, www.arabidopsis.org) and O. sativa genome database (rice.plantbiology.msu.edu).

Identification of TIFY genes

To identify all possible members of the TIFY gene family in Moso bamboo, the Arabidopsis TIFY sequences were first submitted to the Pfam database (http://pfam.sanger.ac.uk) (Finn et al., 2010) to obtain the domain architecture of this family. The TIFY domain, Jas and CCT motifs were represented by Pfam accession numbers PF06200, PF09425 and PF06203, respectively. The corresponding profiles were downloaded and used to search the Moso bamboo Genome Database (www.bamboogdb.org/) by using HMMER (Eddy, 1998). The obtained candidates were further confirmed by the Pfam database. The full length cDNA and EST datasets, comprising of 10,608 low redundant full-length cDNA sequences and approximately 38,000 ESTs from leaf, shoot, and seedling libraries, were employed to obtain evidences for the obtained genes.

Multiple alignments and phylogenetic analysis

The amino acid sequences of TIFY genes were aligned by using MUSCLE (Edgar, 2004) with default parameters, and the phylogenetic tree was then constructed by using MEGA5.2 software with Neighbor-Joining method (Tamura et al., 2011). The Jones-Taylor-Thornton (JTT) model and a discrete Gamma distribution (+G) with 5 rate categories were chosen based on model test. The bootstrap test was carried out with 1,000 iterations.

Analyses of exon-intron structures and protein motifs

The coordinates of exon and intron of TIFY genes of P. edulis were extracted from their corresponding scaffolds and exon-intron structures were illustrated using Gene Structure Display Server (GSDS, http://gsds.cbi.pku.edu.cn/) (Guo et al., 2007). The consensus amino acid sequences of TIFY domain and Jas motif were illustrated by using WebLogo (Crooks et al., 2004).

Evolutionary analyses of the paralogs and orthogos in P. edulis, O. sativa and Brachypodium

Reciprocal BLASTP was carried out to establish orthologous relationship among P. edulis, O. sativa and Brachypodium. The hits threshold values were set as E-value <1e–10, score >200, and positive >70%. The paralogous relationship within P. edulis was also analyzed with more stringent parameters of E-value <1e–50, score >200, and positive >80%. To further confirm the paralogous and orthologous relationship, the nucleotide sequences of coding sequences were also aligned and phylogenetic tree was constructed by using p-distance model. The bootstrap test was carried out with 1,000 iterations.

The synonymous (Ks) and non-synonymous (Ka) substitution rates of paralogous and orthologous were analyzed by Ka_Ks calculator 2.0 (Zhang et al., 2006). Time (million years ago, Mya) of duplication and divergence was calculated using a synonymous mutation rate of one substitutions per synonymous site per year as T = Ks/2λ (λ = 6.5 × 10⁻⁹) (Lynch & Conery, 2000; Yang et al., 2008).

Stress treatment and RNA-seq

To evaluate expression patterns of PeTIFY genes under abiotic stress, dehydration and cold treatments were conducted. The P. edulis plants used in this study were manually planted and grown under natural condition at Linyanshan Experimental Base (N31°0′33.20″, E103°36′51.95″) of Sichuan Agricultural University, Dujiangyan, Sichuan, China. The approximately 20 cM branch containing full young unexpanded leaves of 4–5 centimeters long were cut from five two-year old P. edulis plants in similar growth status. All these samples were collected in the morning at 10, October, 2015 (approximately at 10 o’clock and it was cloudy with temperature of ∼22 °C). For dehydration treatment, the branches were placed on the dry filter paper and treated under room temperature (20 °C and 50% humidity). For cold treatment, the branches were put into a chamber set to 0 °C without light. At 2h and 8h after each treatment, ten individual unexpanded leaves were detached from base and immediately frozen in liquid nitrogen, from which the total RNA was extracted according to the manual of the TRIZOL RNA Kit (TIANGEN, Beijing, China). The same amount of untreated leaves were also sampled and immediately frozen in liquid nitrogen, which were used as control. The qualities and quantities of extracted nucleotide were measured by NanoDrop 2000 Spectrophotometer (Thermo Fisher, USA) and Agilent 2100RNA 6000 Nano kit.

The cDNA library construction and sequencing on Illumina HiSeq™ 4000 platform were performed by Onmath Co. (Chengdu, China), following the manufacturer’s standard protocol. The sequences of 150 bp in length were generated as raw data by pair-end sequencing strategy. The filtered clean reads were mapped to all obtained PeTIFY genes by using TopHat v2.0.9. HTSeq v0.6.1 was used to count the reads numbers mapped to each gene. And then RPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene. RPKM, Reads Per Kilobase of exon model per Million mapped reads, considers the effect of sequencing depth and gene length for the reads count at the same time, and is currently the most commonly used method for estimating gene expression levels (Mortazavi et al., 2008).

Quantitative real-time PCR (qPCR)

The cDNA was synthesized with Oligo18dT primer by using single-stranded cDNA Synthesis Kit (TaKaRa, Dalian, China) following the manufacture introduction. qPCR reactions were conducted employing an BioRad CFX system using the iQ SYBR Green supermix kit (BioRad). The PCR procedure was: pre-incubation at 95 °C for 5 min followed by 40 cycles of denaturation at 95 °C for 15 s, annealing at 60 °C for 15 s, and extension at 72 °C for 15 s. The elongation factor 1a (EF1) (Long et al., 2010; Liang et al., 2012) was used as internal control to quantify the relative transcript level. The primers with a single peak on the melting curve were chosen for qPCR analyses (Table S1). The target and internal control were amplified in separate wells in triplicate. The Cq values were determined automatically by BioRad CFX manager 2.1 (BioRad) and the mean Cq of triplicate was used to calculate the relative level of gene expression by using the 2^−ΔΔct method (Livaka & Schmittgen, 2001). The samples from three independent biological replicates were analyzed and compared to those obtained by RNA-seq.

Identification of TIFY genes in P. edulis

To identify putative TIFY proteins from P. edulis genome, the hidden Markov model (HMM) profiles of the TIFY domain, Jas and CCT motifs were downloaded from Pfam database and utilized to screen the P. edulis genome data by using HMMER. A total of 24 proteins containing a TIFY domain were identified, in which five proteins contain both a TIFY domain and a CCT motif, 17 proteins contain both a TIFY domain and a Jas motif (also known as CCT_2 motif), respectively. Two proteins, PH01000549G0400 and PH01000361G0580, only contain a TIFY domain. The domain or motif of these proteins were re-analyzed by Pfam and showed consistent results with those of HMMER. Additionally, a C2C2-GATA zinc-finger domain (GATA) was found in three out of five proteins containing both a TIFY domain and a CCT motif (Table 1). Using these protein sequences as queries, we searched the EST and full length cDNA of P. edulis. As a result, 16 TIFY genes obtained evidence from expressed sequences (Table 1).

Interestingly, two GATA domains are present in PH01000114G0660, whereas no GATA domains was found in one protein PH01000878G0620, who also has both the TIFY domain and CCT motif. We also noticed that PH01000361G0580 contains a incomplete TIFY domain with a rather lower E-value (9.70E–05). Therefore, to further validate these equivocal gene models, we retrieved their corresponding scaffolds and re-predicted the gene models by GENESCAN (Burge & Karlin, 1998). As results, PH01000114G0660 was predicted to be same as previously presented. Two novel gene models were harvested: PH01000361G0580 was re-predicted to encode a protein of 209 amino acids rather than previous 83 amino acids, containing both the TIFY domain and Jas motif (Table 1, Table S2, and Fig. S1). PH01000878G0620 was predicted to have a coding sequence (CDS) of 5,157 bp, which encodes a protein of 1718 amino acids, containing a TIFY domain, a Jas motif, as well as a GATA domain. The length of genomic sequence of the new PH01000878G0620 is 22,593 bp, which is much longer than previous 2,253 bp (Tables 1, S2, and Fig. S1). Comparing to two homologous proteins in O. sativa and Brachypodium (LOC_Os06g48534.1 and Bradi1g33980), a peptide comprising of 676 amino acids seems to be inserted between CCT motif and GATA domain of PH01000878G0620, and additionally, it is fused with another protein of 794 amino acids at C-terminal domain (Fig. S2). Protein annotation by Pfam and DELTA-BLAST algorithm on NCBI website indicated that the 676 amino acids insertion consists of an UBN2_3 domain (gag-polypeptide of LTR copia-type), a RVT_2 (reverse transcriptase) peptide, and a RNase_HI_RT_Ty1 (Ty1/Copia family of RNase HI in long-term repeat retroelements) peptide (Table S3, Fig. S2). The C-terminal of PH01000878G0620 could be annotated as Arsenite-resistance protein 2 (ARS2), and comprises a DUF3546 (Domain of unknown function) domain (Table S2 and Fig. S2). These results indicate that a putative retrotransposition event may have occurred leading to putative structural variation of PH01000878G0620.

Classification and gene structure

Based on the multiple alignment of PeTIFY protein sequences, a neighbor-joining tree was constructed (Fig. 1). It showed that the 18 TIFY proteins containing a TIFY domain and a Jas motif were clustered together and they were assigned as JAZ subfamily; Five proteins containing a TIFY domain a CCT motif, and a GATA domain were clustered closely and belong to ZML subfamily; One protein only contained a TIFY domain and was clustered separated from others, which was assigned to TIFY subfamily. Compared with some plant species (Bai et al., 2011), P. edulis shares similar gene member to those of monocots, such as O. sativa, B. distachyon, and Sorghum bicolor, but different from those of dicots species.

By drawing the sequence logos, we analyzed conservation of TIFY domain and Jas motif (Fig. 1C). Besides TIFY motif, Gly-12 (position is according to that present in Fig. 1C), Val-14, Val-16, Phe-17, and Lys-24 are conserved crossing the 24 PeTIFYs. Compared with TIFY domain, the Jas motif has more conserved residues (Fig. 1D). The JAZ proteins play important roles in Jasmonates-mediate responses to biotic and abiotic stresses, as well as some important development progresses (Wasternack, 2007). For example, some JAZ proteins are transcriptional repressors of Jasmonates-responsive transcription factors (such as MYC2 and MYC3). With the accumulation of Jasmonates, JAZ will interact with CORONATINE-INSENSITIVE1 (COI1), and trigger the ubiquitination degradation of JAZ proteins mediated by COI1- E3-type ubiquitin ligase SCF (Skp/Cullin/F-box) complex (SCF^COI1), and the de-repression of transcription factors (Yan et al., 2007; Thines et al., 2007; Melotto et al., 2008). Due to the essential role in this interaction and Jasmonates sensitivity, Jas motif is highly conserved in both diocts and monocots (Fig. S3).

The exon/intron structures among the 24 TIFY genes of P. edulis are divergent. Except for five genes containing an uninterrupted open reading frame, all others showed exon/intron structures. It was also found that some genes with close phylogenetic relationship also share conserved exon/intron structure. We identified four genes pairs, PH01000038G0470/PH01000115G0020, PH01000213G1380/PH01001078G0420, PH01000310G0500/PH01000158G0210, and PH01002950G0020/PH01000008G2960, exhibited same number of exons with highly similar length to each other, respectively (Fig. 1E).

Phylogenetic analysis of TIFY family among P. edulis, O. sativa, and Arabidopsis

To analyze the evolutionary relationships of the TIFYs, a Neighbor-Joining tree was constructed based on the alignment of 24 PeTIFYs, 20 TIFYs of O. sativa, and 18 TIFYs of Arabidopsis (Vanholme et al., 2007). According to the phylogenetic tree, all 62 proteins could be separated to several clusters (Fig. 2). All ZML subfamily members were clustered together. PH01000549G0400, LOC_Os02g49970.1 and AT4G32570.1 were grouped closely, which belong to the TIFY subfamily and were close to ZML subfamily. All proteins of JAZ subfamily formed several groups. JAZ I group is composed of four PeTIFYs, three OsTIFYs, and four AtTIFYs, respectively. JAZ II group is close to JAZ I, consisting of nine PeTIFY proteins. No AtTIFY was clustered in this group. Similar results were found in previous studies (Ye et al., 2009; Zhang, You & Chan, 2015), suggesting that group JAZ II might be specific for monocots. Four PeTIFYs were grouped together with eight AtTIFYs and for OsTIFY proteins, which were assigned to group JAZ III. The repredicted PH01000361G0580 was grouped closely with two AtTIFYs and two OsTIFYs, which were assigned as group JAZ IV.

Among group JAZ III, two PEAPOD (PPD) proteins of Arabidopsis were cluster separately. The PPD protein, containing an N-terminal PPD domain, was firstly identified by map-based cloning in ppd mutants (White, 2006), in which PPD1 has been suggested to coordinate tissue growth, modulate lamina size, and limit the curvature of the leaf blade (White, 2006), while the PPD2 was found to interact with the geminivirus AL2 protein and the coat protein promoter (Lacatus & Sunter, 2009). They have been found in dicot species, such as Populus trichocarpa, Mimulus guttatus, Glycine max, Medicago trunculata, and Vitis vinifera, etc (Bai et al., 2011; Zhang et al., 2012), but had not been identified in monocots.

Duplication, orthologous relationship and divergence rates among PeTIFYs, OsTIFYs, and BdTIFYs

P. edulis, O. sativa and B. distachyon are representatives of Bambusoideae, Ehrhartoideae and Pooideae, the three subfamilies of grass (BEP clade), whose phylogenetic relationship was just resolved recently (Wu & Ge, 2012; Peng et al., 2013). It is interesting to investigate the evolutionary relationship of TIFY family between the three species. By reciprocal BLASTP analyses, 15 OsTIFYs (75%) and 12 BdTIFYs (57%) are found to be orthologous to 20 PeTIFYs (83%), respectively (Table S4).

Whole genome duplications, such as tandem and segmental duplications usually give rise to copy numbers of a gene family and have been considered as important drive force of occurrence of novel biological functions. The reciprocal BLASTP analysis revealed that 12 PeTIFYs are involved in six paralgous gene pairs, in which one pair belongs to ZML subfamily, and the other five belong to JAZ subfamily, respectively. We also conducted phylogenetic analyses of three species to support the paralogous and orthologous relationship identified by BLAST analyses. According to the phylogenetic tree, the paralogues or orthologues were clustered together with high bootstrap values (Fig. 3).

In order to evaluate the timing of intragenomic gene duplication events, as well as divergence of orthologoues, the synonymous substitution rate (K_S) was calculate. The paralogous gene pairs exhibited mean Ks around 0.22 (Fig. 3). Then estimated by universal substitution rate of 6.5 ×10⁻⁹ mutations per site per year, the duplications of PeTIFYs may have occurred around 16.7 million years ago (MYA) (Fig. 4). Peng et al. (2013) reported a putative recent whole genome duplication at 7-12 MYA based on analyzing whole genome duplication events of P. edulis. Wu et al. (2015) found that P. edulis genome underwent a large-scale event around 15 MYA by genome-wide analysis of AP2/ERF transcription factors. Our result is similar to these reports and supported that the expansion of PeTIFYs occurred along with the genome duplication event. Among the P. edulis-O. sativa, P. edulis –B. distachyon, and O. sativa- B. distachyon orthologous gene pairs, the mean Ks are approximately 0.51, 0.51, and ∼0.58, indicating that their divergent time was ∼39 MYA, 39 MYA, and 45 MYA, respectively (Fig. 4A). This is consistent with that inferred by whole genome analyses (Peng et al., 2013).

We also calculated the ratios of non-synonymous (Ka) versus synonymous (Ks) substitution rate (Ka/Ks) for duplicated gene-pairs as well as the orthologues of O. sativa and B. distachyon (Fig. 4B). The Ka/Ks ratio is a measure of the selection pressure to which a gene pair is subjected. Ka/Ks <1 means purifying or negative selection, Ka/Ks = 1 stands for neutral selection, and Ka/Ks >1 indicates positive selection (Lynch & Conery, 2000). The Ka/Ks for paralogous gene pair of PeTIFYs are 0.3–0.91 with mean of 0.5. Those for orthologous genes pairs of PeTIFY-OsTIFY, PeTIFY-BdTIFY, and OsTIFY-BdTIFY are 0.23–0.68 with mean of ∼0.38, 0.26–0.59 with mean of ∼0.38, respectively (Fig. 4B). These results indicated that TIFY family in the three species appear to have undergone extensive purifying selection during evolution.

Abiotic stress-induced expression profiling

Plant hormones, such as abscisic acid (ABA), jasmonic acid (JA), salicylic acid (SA), and ethylene (ET) are small signaling molecules involved in stress response. The TIFY family were initially found to be involved in JA signaling pathway. The further genome-wide analyses indicated that most of TIFY genes are regulated by JA and ABA, but not SA and ET (Zhang et al., 2012; Li et al., 2015). As both phytohormones are responsive to biotic and abiotic stresses, the TIFY genes may have versatile roles in stress responses. In Arabidopsis, genes of TIFY family, especially JAZ subfamily, are induced by wounding, herbivory, and pathogen infection (Chung et al., 2008; Demianski, Chung & Kunkel, 2012). Additionally, several AtTIFY genes in Arabidopsis are also found to be upregulated by ozone (Rao et al., 2000), and high salt concentrations (Jiang & Deyholos, 2006). Since bamboo is one of the fast growing and economically important non-timber forest resources in the world, abiotic stress can have enormous negative impact on its production. Therefore, it is meaningful to excavate PeTIFY genes responsive to abiotic stress, such as dehydration and cold.

We performed transcriptome sequencing to evaluate dynamic expression levels of PeTIFY genes under dehydration and cold stresses (Fig. 5A). Approximately four giga bases high quality data (error rate in base calls is less than 0.1%.) for each sample were generated and used to calculate RPKM values of PeTIFYs. Two biological replicates were performed and showed high degree of correlations (Fig. S4). We further validated transcripts levels by using qPCR. As results, most of the genes exhibited similar expression patterns to those harvested by RNA-seq (Fig. S5), indicating high reliability of RNA-seq data in this study.

Approximately 37% of PeTIFYs are highly expressed (log2 RPKM value > 4) in leaf under normal conditions, 46% showed moderate expression levels (log2 RPKM value between 1 and 4), and the remaining 17% are in very low expression level or not expressed (Figs. 5A and 5B). Among the two repredicted genes, PH01000361G0580 seemed not expressed, whereas PH01000878G0620 showed low expression level. Most of the highly expressed PeTIFY genes are of JAZ subfamily, whereas, all ZML subfamily member showed rather low expression level or no expression. As the JAZ proteins are also involved in regulation of some important developmental progresses, such as seed germination, root growth, and flower development (Wasternack, 2007), we proposed that the JAZ genes of P. edulis expressed under normal conditions may function in development process of P. edulis.

Under dehydration and cold stresses, 11 and 13 TIFY genes are up-regulated (RPKM fold change > 2), two genes are down-regulated, while seven genes are not responsive to the two stresses (Figs. 5A and 5C). All PeTIFY genes up-regulated by dehydration could also be up-regulated by cold. Similar phenomenon was found in different plant species. In grape, the majority of the TIFY genes were responsive to osmotic and cold stress, almost all of genes exhibiting responses to cold were also regulated by drought or salinity (Zhang et al., 2012). In apple, 10 of the 22 apple TIFY genes respond to drought or high salinity, and among six were upregulated by both treatments (Li et al., 2015). These results suggested that these genes may have multiple functions in responding to the different abiotic stresses.

Five PeTIFYs, PH01001078G0280, PH01001852G0020, PH01000213G1410, PH01000038G0470, and PH01000115G0020, show low expression levels under normal condition and are quickly responsive to both dehydration and cold, in which PH01001852G0020, and PH01000038G0470 and PH01000115G0020 are homologous to OsTIFY11a and OsTIFY11c of rice, respectively. These two genes function in salt and dehydration tolerance (Ye et al., 2009; Toda et al., 2013). Additionally, PH01001078G0420 and PH01000213G1380 are also quickly respond to stress (Fig. 5A). They are homologous to OsTIFY11d involved in drought tolerance (Ye et al., 2009). The correspondence between expression patterns and homology to function known TIFY genes suggest that they may have similar functions in stress responses. Further work is necessary to functionally characterize the roles of these stress-induced PeTIFY genes in responding to abiotic stress.