Gene expression profile indicates involvement of NO in Camellia sinensis pollen tube growth at low temperature

Illumina sequencing and de novo assembly

To study the transcriptome expression of pollen tubes, C. sinensis pollen tubes were cultured under the following conditions: 25 °C (CK), 4 °C (LT) and 25 μM DEA NONOate (NO). Three libraries (CsPT-CK, CsPT-LT, and CsPT-NO) were designed for RNA-Seq using an Illumina HiSeq^TM2000 genome analyzer, and all of the experiments were replicated three times. After cleaning and performing quality checks, 55.5, 52.4, and 51.3 million clean reads were obtained for CsPT-CK1, CsPT-CK2 and CsPT-CK3 (Accession number SRR3270364, SRR3270376 and SRR3270829 for library CsPT-CK1, CsPT-CK2 and CsPT-CK3), respectively; 50.2, 51.8, and 51.4 million clean reads were obtained for CsPT-LT1, CsPT-LT2 and CsPT-LT3 (Accession number SRR3270928, SRR3270974 and SRR3270993 for library CsPT-LT1, CsPT-LT2 and CsPT-LT3), respectively; and 51.9, 54.5, and 53.7 million clean reads were obtained for CsPT-NO1, CsPT-NO2 and CsPT-NO3 (Accession number SRR3270997, SRR3271001 and SRR3271002 for library CsPT-NO1, CsPT-NO2 and CsPT-NO3), respectively. The Q20 percentage was > 97.77 % in each library as shown in Table 1. The clean reads were assembled into 60,611, 61,246 and 63,117 contigs with mean lengths of 333, 335 and 331 nt for CsPT-CK, respectively; 56,102, 58,622 and 60,547 contigs with mean lengths of 343, 340 and 333 nt for CsPT-LT, respectively; and 59,750, 58,884 and 59,550 contigs with mean lengths of 331, 335 and 333 nt for CsPT-NO, respectively. These contigs were further assembled into 39,726, 40,440 and 41,626 unigenes with mean lengths of 648, 650 and 660 nt for CsPT-CK, respectively; 36,993, 39,070 and 39,439 unigenes with mean lengths of 650, 664 and 657 nt for CsPT-LT, respectively; and 39,514, 38,298 and 39,061 unigenes with mean lengths of 647, 650 and 639 nt for CsPT-NO, respectively (Table 2). These unigenes size distribution was displayed in Additional file 1: Figure S1. As shown in Fig. 1, the Pearson product moment coefficient was > 99.65 % in each library. All above data indicated that the quality of throughput and sequencing was high enough for further analysis.

Sample	Total raw reads	Total clean reads	Total clean nucleotides (nt)	Q20 percentage	N percentage	GC percentage
CsPT-CK-1	54,971,748	51,933,988	4,674,058,920	97.92 %	0.00 %	45.44 %
CsPT-CK-2	57,928,046	54,556,816	4,910,113,440	97.92 %	0.00 %	45.23 %
CsPT-CK-3	56,727,452	53,686,756	4,831,808,040	97.93 %	0.00 %	45.39 %
CsPT-LT-1	58,915,630	55,506,666	4,995,599,940	97.86 %	0.00 %	45.64 %
CsPT-LT-2	55,468,880	52,434,854	4,719,136,860	97.85 %	0.00 %	45.39 %
CsPT-LT-3	54,793,794	51,330,908	4,619,781,720	97.77 %	0.01 %	45.35 %
CsPT-NO-1	53,553,772	50,257,976	4,523,217,840	97.84 %	0.01 %	45.20 %
CsPT-NO-2	55,260,058	51,851,106	4,666,599,540	97.82 %	0.01 %	45.40 %
CsPT-NO-3	54,754,898	51,364,466	4,622,801,940	97.83 %	0.01 %	45.30 %

CK control, LT 4 °C treatment, NO NO donor DEANONOate treatment. All of the experiments were replicated three times. Q20 percentage shows the percentage of nucleotides with a quality value higher than 20; N percentage represents the percentage of unknown nucleotides in clean reads; and GC percentage indicates the percentage of guanidine and cytosine nucleotides among total nucleotides

	Sample	Total number	Total length(nt)	Mean length(nt)	N50	Total consensus sequences	Distinct clusters	Distinct singletons
Contig	CsPT-CK-1	60,611	20,192,162	333	571
	CsPT-CK-2	61,246	20,491,956	335	581
	CsPT-CK-3	63,117	20,892,307	331	561
	CsPT-LT-1	56,102	19,226,389	343	599
	CsPT-LT-2	58,622	19,960,405	340	595
	CsPT-LT-3	60,547	20,132,613	333	572
	CsPT-NO-1	59,750	19,797,938	331	559
	CsPT-NO-2	58,884	19,734,167	335	575
	CsPT-NO-3	59,550	19,805,943	333	568
Unigene	CsPT-CK-1	39,726	25,741,808	648	1165	39,726	14,095	25,631
	CsPT-CK-2	40,440	26,270,780	650	1152	40,440	14,717	25,723
	CsPT-CK-3	41,626	27,483,327	660	1204	41,626	15,310	26,316
	CsPT-LT-1	36,993	24,057,930	650	1182	36,993	13,290	23,703
	CsPT-LT-2	39,070	25,937,230	664	1203	39,070	14,224	24,846
	CsPT-LT-3	39,439	25,906,326	657	1191	39,439	14,301	25,138
	CsPT-NO-1	39,514	25,545,806	647	1173	39,514	14,211	25,303
	CsPT-NO-2	38,298	24,889,732	650	1169	38,298	13,676	24,622
	CsPT-NO-3	39,061	24,971,317	639	1153	39,061	13,923	25,138
	All	45,432	46,435,079	1022	1574	45,432	22,812	22,620

Total Consensus Sequences represents all of the assembled unigenes. Distinct Clusters represents the clustered unigenes, with the same cluster containing high-similarity (more than 70 %) unigenes that may be from the same gene or homologous genes. Distinct Singletons represents unigenes from a single gene

Annotation of predicted proteins

We annotated 36,097 unique sequences on the basis of a BLAST search of six public databases: NCBI non-redundant (NR) database, mysql-nt (NT) database, Swiss-Prot protein database, Kyoto Encyclopedia of Genes and Genomes (KEGG) database, Cluster of Orthologous Groups of proteins (COG) database and Gene Ontology (GO) database. As shown in Fig. 2, we annotated 33,716 unique sequences using the NR database as a reference. Based on these NR annotations, 48.4 % of the sequences had a very strong homology to available plant sequences (E-value < 10⁻⁶⁰), 20.4 % showed a strong homology (10⁻⁶⁰ < E-value < 10⁻³⁰), and an additional 31.3 % showed homology (10⁻³⁰ < E-value < 10⁻⁵) (Fig. 2a). Similarity distributions were comparable, with 31.5 % of the sequences presenting similarities higher than 80 % and 69.5 % presenting similarities of 17–80 % (Fig. 2b). Furthermore, we found 43.7 % of the sequences matching homologues from Vitis vinifera, with additional hits to Amygdalus persica (10.1 %), Lycopersicon esculentum (10.1 %), Ricinus communis (7.9 %), and Populus balsamiferas ubsp. trichocarpa (6.7 %) (Fig. 2c).

Functional classification by COG

We searched all unigenes against the COG database for functional predictions and classifications. Our query assigned 13,029 of the 33,716 sequences that presented NR hits to COG classifications (Fig. 3). The COG-annotated putative proteins were functionally classified into at least 25 molecular families, including cellular structure, biochemistry metabolism, molecular processing, and signal transduction. Among them, the cluster “General function prediction” represented the largest group (4100; 16.96 %), followed by “Transcription” (2073; 8.57 %), “Replication, recombination and repair” (2064; 8.54 %), “Posttranslational modification, protein turnover and chaperones” (1992; 8.24 %), “Signal transduction mechanisms” (1713; 7.09 %), “Translation, ribosomal structure and biogenesis” (1560; 6.45 %), and “Carbohydrate transport and metabolism” (1503; 6.22 %); Only several unigenes were assigned to “Extracellular structure”.

DEGs and GO classification

We filtered the DEGs (Additional file 2: Table S1, Additional file 3: Table S2 and Additional file 4: Table S3) using the absolute values of log₂Ratio > 1 and probability > 0.7 as the thresholds according to Zhao et al. [31]. Comparing the library CsPT-CK with CsPT-LT (CK-VS-LT) showed that there were 766 DEGs (243 up-regulated and 523 down-regulated genes, 243/523). Moreover, we found 497 (222/275) and 929 (545/384) DEGs in the comparison between the library CsPT-CK and CsPT-NO (CK-VS-NO) and between the library CsPT-LT and CsPT-NO (LT-VS-NO), respectively (Fig. 4a). Among these DEGs, 357, 247 and 496 genes were specifically expressed in the CK-VS-LT, CK-VS-NO and LT-VS-NO, respectively. In addition, we observed 101, 284 and 125 genes co-expressed in CK-VS-LT and CK-VS-NO, CK-VS-LT and LT-VS-NO and CK-VS-NO and LT-VS-NO, respectively, and 24 genes were co-expressed in the three comparisons (Fig. 4b).

GO classifies genes on the basis of three functional ontologies, that is, Molecular function, Cellular component and Biological process. We used Blast2GO [32] to obtain GO annotations and Web Gene Ontology Annotation Plot (WEGO) [33] to perform GO functional classifications. For the CK-VS-LT, CK-VS-NO and LT-VS-NO libraries, 262 out of 766 DEGs (262/766), 156 out of 497 DEGs (156/497) and 328 out of 929 DEGs (328/929) could be assigned a GO classification, respectively (Additional file 5: Table S4, Additional file 6: Table S5 and Additional file 7: Table S6). The 262, 156 and 328 DEGs were assigned to 1957, 1239 and 2474 GO classifications in the CK-VS-LT, CK-VS-NO and LT-VS-NO libraries, of which 952, 593 and 1219 as “Biological process”, 705, 459 and 920 as “Cellular component” and 300, 187 and335 as “Molecular function”, respectively (Fig. 5). As shown in Fig. 5, the majority of DEGs in the “Biological process” category were associated with cellular processes, metabolic processes and single-organism processes, most DEGs in the “Cellular component” category were associated with cells, cell parts and organelles, and most in the “Molecular function” category were associated with binding and catalytic activity. Furthermore, between the CK-VS-LT and CK-VS-NO libraries, certain CK-VS-LT DEGs annotated with the “Biological process” category were related to biological adhesion and locomotion, whereas several CK-VS-NO DEGs annotated with the “Molecular function” category were associated with enzyme regulator activity and protein binding transcription factor activity.

Cluster analysis of DEGs co-expressed in CK-VS-LT and CK-VS-NO

Genes are usually functionally correlated to those with similar expression patterns. In order to understand the expression patterns of the 125 co-expressed DEGs (Additional file 8: Table S7) in the CK-VS-LT and CK-VS-NO libraries, a cluster analysis of the gene expression patterns between these two libraries were performed. We classified these genes into three groups (Fig. 6). The largest group (Group 1) contained 71 genes (56.8 %) (from Unigene 11102_All to Unigene 20419_All) (Additional file 8: Table S7) that were up-regulated in both CK-VS-LT and CK-VS-NO. This group mainly contained genes encoding proteins related to oxidation-reduction reactions and amino acid metabolism, such as cytochrome c oxidase, arginine/serine-rich protein and serine/threonine-protein kinase. The second largest group (Group 3) included 52 genes (41.6 %) (from Unigene 10170_All to Unigene 8447_All) that were down-regulated in both CK-VS-LT and CK-VS-NO. This group mainly contained genes encoding proteins related to plant hormone metabolism and signaling pathways and TFs, such as ethylene receptor, ABA 8’-hydroxylase (ABA8ox), MADS-box protein and CaM-binding transcription activator (CAMTA). Finally, two genes (Group 2) were down-regulated in CK-VS-LT but up-regulated in CK-VS-NO, and one of them (Unigene 11403_All) was annotated as a U-box domain-containing protein.

DEGs identified from pollen tubes treated with low temperature and NO

Quantitative real-time PCR (qRT-PCR) validation of DEGs in CK-VS-LT and CK-VS-NO

To confirm the gene expression patterns shown in the sequencing data, we performed a qRT-PCR analysis on nine randomly selected genes in pollen tubes exposed to low-temperature stress, 25 μM DEA NONOate (NO donor) or to low temperature with 200 μM 2-(4-carboxyphenyl)-4,4,5,5-tetramethylimidazoline-1-oxyl-3-oxide (cPTIO) (NO scavenger). According to our preliminary analysis and previous reports, we speculated that Ca²⁺ signaling and TFs played an important role in this process [11, 25]. Therefore, we selected CAMTA (Unigenes 2215), MYB (Unigenes 17702), ZFP (Unigenes 21920), AP2 (CL2776.Contig1) and PP2C (CL580.Contig48) genes. Furthermore, Phospholipase D (PLD) accumulated highly when Chorispora bungeana was exposed to low temperature [34] and PLD activation in C. bungeana exposed to low-temperature stress was related to a drop in Ca²⁺ content in membrane fractions [35]. Moreover, TOPLESS (TPL) was involved in mutual adjustment between hormones via NO in Petunia hybrida [36]. Additionally, genes encoding a non-specific serine/threonine protein kinase were down-regulated expressed during response to low temperature in two grapevine varieties [37]. Recent research indicated that phytochrome B (PHYB) and PHYA antagonistically regulated tomato low-temperature tolerance, in a process involving far-red light-induced PHYA activation to induce ABA signaling and eventually lead to cold response [38]. Here, we also randomly selected PLD (Unigenes 22243), TPL (Unigenes 8699) and serine/threonine protein kinase (CL10.Contig2) in CK-VS-LT and CK-VS-NO, and the PHYB (CL2425.Contig1) in CK-VS-LT for our qRT-PCR assay. The expression of Unigene 8699, CL10.Contig2 and CL580.Contig48 was up-regulated, whereas Unigenes 2215, 17702, 21920 and 22243 and CL2425.Contig1 and CL2776.Contig1 were down-regulated expressed in CK-VS-LT and CK-VS-NO (Fig. 7). These analyses support our DEGs data, and the high confirmation rate further confirmed the reliability of the data. In addition, application of 200 μM cPTIO under low-temperature stress rescued the changes in the expression levels of the DEGs co-expressed in the three treatments, such as Unigenes 2215, 8699, 17702, 21920 and 22243; CL10.Contig2 and CL580.Contig48. These results indicated that the seven genes intervened in pollen tube response to low temperature, possibly via the regulation of NO signaling.

Comparison of unigenes between C. sinensis leaves and pollen tubes

According to previous reports, Illumina sequencing of leaves transcriptome from C. sinensis cv. ‘Longjing43’ revealed an average transcript size of 356 bp with a N50 of 529 bp [5]. In this study, we reported that the average transcript size was > 331 bp, and the N50 was > 559 bp in every transcriptome of pollen tubes from C. sinensis cv. ‘Longjingchangye’. In addition, the previous study showed the annotation of 53,201 unigenes based on BLASTx searches against NR (51,883 unigenes), UniRef90 (52,217 unigenes), the Arabidopsis Information Resource (TAIR, version 10) (42,969 unigenes), KEGG (15,086 unigenes) and COG (29,356 unigenes) databases in tea plant leaves [5]. Among 71,289 transcripts based on NR annotation, 30,115 transcripts matched 80,176 GO term annotations [5]. Most of the unigenes both from leaves and pollen tubes are categorized as the cell and cell part (cellular component), binding and catalytic activity (molecular function) and cellular process and metabolic process (biological process) (Additional file 15: Figure S2). Interestingly, envelope (cellular component), auxiliary transport protein, transcription regulator (molecular function) and anatomical structure formation, cell killing, death, pigmentation, viral reproduction (biological process) were only annotated in GO categories of tea leaves transcriptome (Additional file 15: Figure S2A). However, cell junction, extracellular matrix, extracellular matrix part, membrane, membrane part, nucleoid, symplast (cellular component), channel regulator, metallochaperone, nucleic acid binding transcription factor, nutrient reservoir, protein binding transcription factor, protein tag, receptor (molecular function) and regulation of biological (biological process) were only annotated in GO categories of pollen tube transcriptome (Additional file 15: Figure S2B). These differences between transcriptomes from leaves and pollen tubes might result from the diversity of reproductive and vegetative tissues of plant [6], further confirming that reproductive tissue responds to low temperature stress using distinct mechanisms from vegetative tissues [7].

Plant hormone metabolism and signaling pathways

Plant hormones play essential roles in response to abiotic stress through their involvement in complicated signaling networks and physiological processes. For example, ABA has been shown to accumulate as a protective factor against cold injury in maize [39], and previous investigations have demonstrated that NO enhances tolerance to drought stress via promoting ABA-induced stomata closure [40, 41]. In our study, we identified one ABA8ox gene (Unigene14044_All) and one PP2C gene (CL580.Contig48_All) in the CK-VS-LT (Table 3) and CK-VS-NO libraries (Table 4). As shown in previous reports, NO-induced ABA content reduction in Arabidopsis correlates with the modulation of ABA8ox protein expression [42]. Additionally, ABA signaling is negatively regulated by NO in germination and early seedling growth in Arabidopsis through the S-nitrosylation of sucrose non-fermenting 1 (SNF1)-related protein kinase 2.2 (SnRK2.2) and SnRK2.3 [43]. As a negative regulator, PP2C inactivates SnRK2, thereby regulating ABA signals [44]. These reports and our studies indicate that NO might participate in ABA catabolism and signaling to regulate pollen tube elongation under cold stress through the genes ABA8ox (Unigene14044_All) and PP2C (CL580.Contig48_All). In addition, the plant steroidal hormone BRs has been found to promote growth and cell elongation and protect plants against abiotic stresses. Reports have shown that BRs exert positive effects on the mitigation of oxidative damage by enhancing the antioxidant defense system in Chorispora bungeana under cold stress [45]. BRs are recognized through an active complex of a leucine-rich repeat (LRR) receptor-like kinase (RLK) brassinosteroid-insensitive 1 (BRI1) and BRI1-associated kinase 1 (BAK1). BRs bind to the extracellular LRR domain of BRI1 and induce a phosphorylation-mediated cascade to modulate gene expression [46]. Meanwhile, BAK1 is involved in the heterodimerization and endocytosis associating with BRI1 to promote BRs-induced signal transducing instead of involved in binding to BRs. [47]. We identified two BRI1 genes (Unigene9712_All and Unigene4038_All) and one BAK1 gene (CL4145.Contig1_All) in CK-VS-LT (Table 3) as well as one BRI1 gene (Unigene4038_All) in CK-VS-NO (Table 4). Current studies have revealed that NO concentration in the root cells is increased by BRs, which is required for BR-induced alterations of the root architecture in Arabidopsis [48]. This finding suggests that BRs are involved in pollen tube resistance to low temperature through Unigene9712_All, Unigene4038_All and CL4145.Contig1_All and indicate that NO might participate in this process through Unigene4038_All. Moreover, BRs could overcome the inhibition of seed germination induced by ABA and increase the tolerance of Arabidopsis seedlings to cold [49]. BR-induced NO production and NO-activated ABA biosynthesis are significant mechanisms for BR-increased tolerance to water-stress [50]. Combining these reports and our findings, BRs and ABA might interact in NO signal transduction during the pollen tube responding to low-temperature stress, although the underlying mechanisms require further research.

Ethylene (ET) is a gaseous plant hormone, and its signaling pathways participate in plant responses to low temperature [51, 52]. Previous reports indicated that ET signaling negatively regulates Arabidopsis tolerance to cold stress by inhibiting the expression of CBF and type-A ARR genes [53]. Ethylene response sensor 1 (ERS1) has been identified as one of the membrane-located receptor proteins in the ET signaling pathway [54]. The level of ET is usually low in normal conditions, and its response is suppressed by ERS1 via activating the downstream negative regulator CTR1 (Constitutive Triple Response 1). Furthermore, EIN2 (ethylene insensitive 2) protein acting as a critical positive regulator of ET signaling has been found to be inhibited by CTR1 in Arabidopsis [55]. As downstream of EIN2, the plant-specific transcription factors EIN3 and EIN3-like1 (EIL1) are sufficient and necessary to stimulate ethylene-response genes expression and to modulate the ethylene-related responses of plants [56]. EIN3-Binding F-box-1 (EBF1) and EBF2 tightly regulate the stability of EIN3 at the protein level [57], and EBF2 transcription is in turn activated by EIN3, thereby forming a negative feedback loop [58]. Interestingly, recent work has shown that ERS1 and EBF1/2 are differentially expressed in Papaya fruit ripening disorders caused by cold injury [59]. In our study, two ERS1 genes (Unigene4486_All and Unigene22534_All) and two EBF1/2 genes (CL2610.Contig1_All and Unigene10222_All) were down-regulated between CK and LT (CK-VS-LT) (Table 3). These findings imply that ERS1 and EBF1/2 are involved in pollen tube resistance to low temperature. Although few reports on the interaction between the two gaseous hormones NO and ET during low-temperature stress are available, studies have suggested that NO and ET have a complicated relationship in plants during water logging stress [60], salt stress [61], UV-B stress [62] and hypersensitive responses [63]. One ERS1 gene (Unigene22534_All) was also identified in the CK-VS-NO library (Table 4). It is tempting to speculate that ERS1 (Unigene22534_All) responds to low temperature via the NO signaling pathway and the other ERS1 genes in C. sinensis pollen tubes (Unigene4486_All and EBF1/2 (CL2610.Contig1_All and Unigene10222_All) respond to low temperature in a NO-independent manner.

Auxin is an essential morphogenetic signal participating in the regulation of cell identity during all the developmental processes of plants, and auxin signaling pathways constitute a critical component of mechanisms for plant tolerance to abiotic stresses. Substantial evidence has demonstrated a tight link between auxin responses and cold stress [64, 65]. For example, cold stress affects auxin polar transport via the selective inhibition of the intracellular trafficking of numerous proteins in Arabidopsis, including the auxin efflux carriers [66]. Two auxin-repressed protein genes (CL2941.Contig1_All and CL2941.Contig2_All), one auxin response factor gene (Unigene7793_All), one indole-3-acetic acid (IAA)-amido synthetase gene (Unigene13496_All) and one IAA-induced protein gene (Unigene21651_All) were identified in the CK-VS-LT library in this study (Table 3). Interestingly, as shown in a recent report, auxin participates in regulating pollen tube tip growth [67]. Therefore, it is reasonable to speculate that auxin intervenes in pollen tube responding to low-temperature stress in C. sinensis. Furthermore, Markus et al. [68] reported that the release of IAA and NO from peroxisomes is a spatially and temporally coordinated process in root formation in maize and Arabidopsis, thereby indicating that auxin and NO interact closely in plant developmental processes. We identified one IAA-amido synthetase gene (Unigene17743_All) and one IAA-induced protein gene (Unigene21651_All) in CK-VS-NO (Table 4). Unigene21651_All was simultaneously detected in both CK-VS-LT and CK-VS-NO libraries, indicating that IAA-induced protein (Unigene21651_All) might be related to the NO-mediated pathway in C. sinensis pollen tube growth in responding to low-temperature stress.

Additionally, gibberellins (GAs) are another class of plant hormones that play critical roles in modulating plant growth and development in response to environmental cues. Exposure of Arabidopsis seedlings to cold stress triggers a reduction in bioactive GA, which restrains root growth [69]. Previous studies have shown that the nuclear-localized growth-repressing DELLA proteins are central components of the GA-signaling pathway. These proteins accumulate through a posttranslational mechanism mediated by the decrease in bioactive GA via up-regulating the transcript level of GA 2-oxidase (GA2ox) gene, and this in turn increases the transcript levels of GA20ox and GA3ox gene by a feedback mechanism [70, 71]. The binding of bioactive GA to GA insensitive dwarf 1 (GID1), a pocket of the GA receptor, promotes the interaction between GID1 and DELLAs and forms the GA–GID1–DELLA complex [72]. This model has been employed to reveal the molecular mechanisms for GA perception. In this study, two DELLA genes (Unigene15959_All and Unigene10398_All), one GA2ox gene (Unigene11055_All), one GA20ox gene (Unigene1092_All) and one GID1 gene (Unigene9129_All) showed a significantly different expression in CK-VS-LT (Table 3). These findings suggest that GA participates in pollen tube response to low-temperature stress, whereas none of the GA signals or synthetic metabolism genes was detected in the DEGs of CK-VS-NO. In addition, María et al. [73] speculated that high levels of NO content could influence the increase of DELLA activity through the disappearance of PIN-FORMED 1 (PIN1) and subsequently inhibit cell elongation in the elongation-differentiation zone in Arabidopsis. However, the PIN1 gene is not detected in CK-VS-NO. Therefore, we hypothesize that NO might regulate C. sinensis pollen tube growth during responses to low temperature stress through a GA-independent signaling network.

Transcription factors (TFs)

Plants have evolved complex and highly efficient regulatory networks to respond and adapt to cold stress, and TFs play crucial roles in these regulatory processes. For example, the rice R2R3-type MYB gene, OsMYB2, plays a positive regulatory role in low-temperature stress tolerance [74]. Similarly, the soybean typical Cys2/His2-type (C2H2-type) zinc finger GmZF1 enhances cold stress tolerance by regulating gene expression in transgenic Arabidopsis [75]. In our study, two MYB genes and three ZFP genes were identified in CK-VS-LT as shown in Table 3. The results reveal that MYB and ZFP are involved in pollen tube elongation in response to low-temperature stress. Besides, NO action inhibits the DNA-binding of AtMYB2 through the S-nitrosylation of Cys53 [76]. Furthermore, the thiols of zinc-sulfur clusters can be S-nitrosated by NO, which reversibly disrupts zinc finger structures and presents a molecular mechanism for regulating ZFP TFs [77]. One MYB gene (Unigene17702_All) and one ZFP gene (Unigene21920_All) were detected among the comparisons of the CK-VS-LT and CK-VS-NO libraries (Tables 3 and 4), implying that MYB (Unigene17702_All) and ZFP (Unigene21920_All) are related to the NO signaling pathway during low-temperature stress responses in C. sinensis pollen tubes. Moreover, NO negatively regulates the DNA binding of MYB2 by S-nitrosylation in the ABA-dependent cold signaling pathway [25]. Combining this information with the results for ABA and NO, we speculate that MYB (Unigene17702_All) might participate in C. sinensis pollen tube growth during low-temperature stress via an ABA-dependent pathway that acts as a NO signaling pathway downstream regulator.

Pollen tube elongation is a tip growing process that can be affected by TFs. For example, the pollen-specific MIKC* class of MADS-box TFs is necessary for pollen maturation and tube elongation and controls a transcriptional switch that directs pollen maturation in Arabidopsis [78]. Furthermore, a large number of reports reveal that calmodulin (CaM) is a key element in pollen tube elongation and orientation and dependent on a complicated crosstalk process between multiple pathways [79–81]. The CaM-interacting protein group CAMTA has been identified in the plant defense signaling cascade [82]. Moreover, MADS-box [83] and CAMTA [84] TFs take part in the regulation of plant responses to low temperature. We detected five MADS-box genes and one CAMTA gene in CK-VS-LT (Table 3). These studies and our results indicate that CAMTA and MADS-box TFs participate in pollen tube elongation under cold stress. In CK-VS-NO, one MADS gene and two CAMTA genes were also detected, and in CK-VS-LT and CK-VS-NO, one MADS gene (CL1553.Contig1_All) and one CAMTA gene (Unigene2215_All) were differentially expressed (Tables 3 and 4). Previous experiments have demonstrated that NO modulates Ca²⁺ signaling in lily and Arabidopsis pollen tubes through directly imaging the concentration of cytosolic free Ca²⁺ during NO-induced tube re-orientation [26, 27]. We speculate that the CAMTA gene (Unigene2215_All) might participate in the regulation of Ca²⁺ by NO during pollen tube elongation under cold stress. Additionally, although directly relevant reports are not available on the crosstalk between MADS-box and NO during pollen tube growth, arginine is the substrate of NO synthases [85], and the metabolism of arginine is regulated by MADS-box [86]. Thus, it is tempting to suggest that the MADS gene (CL1553.Contig1_All) might influence C. sinensis pollen tube elongation under low temperature via the NO synthesis pathway.

In addition to the above TFs, we also identified two ERF genes, one WRKY gene and two bHLH genes in CK-VS-LT (Table 3) as well as one ERF gene and two bHLH genes in CK-VS-NO (Table 4). Using forward and reverse genetics, the ERF [51], WRKY [87] and bHLH [88] TFs were isolated and identified as the molecular switches of cold-responsive signaling networks in plants. Although these three TF family genes were found in CK-VS-LT and CK-VS-NO, none were simultaneously found in both CK-VS-LT and CK-VS-NO libraries. Thus, the complex networks and mechanisms of NO, ERF, WRKY, bHLH and low-temperature stress require further research.

Vesicle polarized trafficking and cell wall biosynthesis related genes

Pollen tube growth is a highly polarized process which is dependent on the polarized trafficking of vesicles containing cell wall materials to the tip region to establish a restricted growth zone [89]. Previous reports have shown that secretory vesicle formation in root hairs is impaired in the double phosphatidylinositol-4-kinase (PI4K) Arabidopsis mutant pi4kβ1/β2 because of the disruptions to phosphatidylinositol-4-phosphate (PI4P) gradient and polarized vesicle trafficking [90]. Furthermore, PI 4-phosphate 5-kinase (PIP5K) isoforms can catalyze the PI 4, 5-bisphosphate (PIP2) synthesis, which eventually regulates vesicle trafficking [91]. Additionally, Feng et al. [92] reported that vesicle fusion and recycling in the tip region of pollen tube apex might be mediated by vesicle-associated membrane protein 726 (VAMP726) in Petunia inflata. Reports have indicated that the vesicle trafficking process is accelerated in response to low temperature stress in Arabidopsis [93]. Here, we identified four PI4K genes, one PIP5K gene and two VAMP genes in CK-VS-LT (Additional file 9: Table S8), implying that these genes modulate pollen tube elongation during the response to cold stress. We also found that one of the PI4K genes (Unigene13861_All) is also differentially expressed in CK-VS-NO (Additional file 10: Table S9). Interestingly, NO involvement in pollen tube growth regulation has been reported to occur via the modulation of vesicle polarized trafficking [29], and NO is also required for vesicle formation and trafficking in root hair tip growth in Arabidopsis [94]. These results indicate that NO may modulate C. sinensis pollen tube elongation under low-temperature stress by mediating vesicle formation and trafficking via the PI4K gene (Unigene13861_All). Additionally, cell wall construction is among the most important mechanisms in pollen tube tip growth [95] and could be regulated by NO in Pinus bungeana [29]. COBRA-like genes are important factors for secondary cell wall development in the root growth zone of Arabidopsis [96]. We also identified one COBRA-like gene (Unigene10170_All) in the DEGs between CK-VS-LT and CK-VS-NO (Additional file 9: Table S8 and Additional file 10: Table S9). This COBRA-like gene (Unigene10170_All) might participate in cell wall construction in C. sinensis pollen tubes during low temperature stress through the NO signaling pathway.

Ubiquitination machinery of the ubiquitin system

The protein ubiquitination pathway is one of the most effective processes for performing versatile post-translational modifications that mediate the growth and development of all eukaryotic species. Ubiquitylated proteins are targeted for degradation by the 26S proteasome which is a proteolytic machine that degrades the target and allows ubiquitin moieties to be reused. Ubiquitin–26S proteasome system (UPS) widely intervenes in plant growth and development during tolerance to stresses by modulating proteins activity or affecting their localization or stability [97]. The UPS depends on the activity of ubiquitin activating enzymes (E1s), ubiquitin conjugating enzymes (E2s) and ubiquitin ligases (E3s). E3s are recognized as the critical factors that define substrate specificity [98] and are divided into four major types: Really Interesting New Gene (RING)-type, cullin–RING ligases (CRLs), Homology to E6-Associated Carboxyl-Terminus (HECT)-type and U-box-type. Furthermore, U-box-type E3s transfer ubiquitin tags directly from E2s-Ubs to target proteins [99]. Additionally, CRLs can consist of a substrate-recruiting protein, such as F-box [100] and broad complex, tramtrack, bric-a-brac/pox virus and zinc finger (BTB/POZ) [101]. In our study, we identified one E2 gene, four E3 genes, one 26S proteasome gene, four U-box genes, nine F-box genes and three BTB/POZ genes in the CK-VS-LT library (Additional file 11: Table S10). UPS modulates the levels of stress hormone biosynthesis and secondary messengers as well as the abundance of regulatory proteins that may be accumulated resulting from exposure to abiotic stress [102], which implies that the UPS plays a fundamental role in C. sinensis pollen tube elongation under low temperature stress. Moreover, most E3s have been linked to abiotic stress responses through their modulation of various stress signaling processes [98]. Previous reports have shown that multiple E3s regulate ABA-dependent stress signaling [103]. Combined with our findings on ABA and cold stress, this information indicates that the UPS is involved in an ABA-dependent signaling pathway and might influence other hormones and TFs through ABA-dependent signaling. Additionally, certain genes encoding E3 proteins are involved in NO-dependent modulation in response to cryptogein in Arabidopsis [104], indicating that NO and UPS interact closely in plants. In CK-VS-NO, four E3 genes, one 26S proteasome gene, two U-box genes, one F-box gene and one BTB/POZ protein gene were identified (Additional file 12: Table S11). Additionally, we also identified one E3s (Unigene6536_All), one 26S proteasome (Unigene21945_All), one U-box (Unigene11403_All) and one BTB/POZ (Unigene17396_All) gene co-expressed in the DEGs of CK-VS-LT and CK-VS-NO. These results indicate that the E3s (Unigene6536_All), 26S proteasome (Unigene21945_All), U-box (Unigene11403_All) and BTB/POZ (Unigene17396_All) might participate in an NO-dependent signaling pathway in C. sinensis pollen tubes under low temperature stress through their influence on the modulation of the UPS, which alters the regulatory proteins levels.

Species-specific secondary metabolites pathways

Major secondary metabolites are fundamental components that define the flavor, taste and health benefits of tea, such as flavonoids, caffeine and theanine. Flavonoids contain flavones, flavonols, isoflavones, flavanones, flavanols, and anthocyanidins, which have various functions in plant tolerance to stresses, including defense against phytopathogens, protection from ultraviolet light and antioxidant activity [105]. Furthermore, the accumulation of catechins in C. sinensis is very susceptible to low temperature [106]. In addition, caffeine is a purine alkaloid widely used as stimulant and pharmaceutical component, and mainly synthesized in young leaves in C. sinensis via a typical biosynthetic pathway that includes steps of purine synthesis and modification [107]. Theanine is emerged as a unique free amino acid that confers the unique “umami” taste of tea, and it occupies about half of the total free amino acids in tea [108]. It has been revealed that theanine biosynthesis starts from glutamine and pyruvate and includes further synthetic steps in buds, leaves and roots from C. sinensis [109, 110]. In the present study, some genes related to flavonoids, caffeine and theanine biosynthesis pathways were down-regulated in CK-VS-LT (Additional file 13: Table S12). In addition, we also identified some genes in CK-VS-NO (Additional file 14: Table S13). And one caffeine synthase 1 (Unigene1362_All) and glutamate receptor gene (CL4694.Contig1_All) were co-expressed in the comparisons of CK-VS-LT and CK-VS-NO. These results reveal that caffeine and theanine mediate cold stress tolerance via NO signal pathway by caffeine synthase 1 (Unigene1362_All) and glutamate receptor (CL4694.Contig1_All).

Plant materials

Mature pollen was collected from C. sinensis (L.) O. Kuntze cv. ‘Longjingchangye’ growing in Dr. Sun Yat-sen’s Mausoleum Tea Factory, Nanjing, China (http://www.yhc1958.com/index.html). Then the pollen was pre-cultured in a control liquid medium (5 % sucrose, 0.05 % Ca(NO₃)₂  · 4H₂O, 0.01 % H₃BO₃, 5 % PEG 4000 and 30 mM MES, pH 6.0) at 25 °C in darkness for 30 min [111]. The pre-cultured samples were incubated at 25 °C for 1 h as control (CK treatment). The pre-cultured pollen was incubated at 4 °C for 1 h as low-temperature treatment (LT treatment). In addition, adding 25 μM NO donor DEA NONOate to the pre-cultured pollen and incubated it at 25 °C for 1 h as NO treatment (NO treatment). For NO scavenger treatment, the pre-cultured pollen was incubated in the medium with 200 μM cPTIO at 4 °C for 1 h (LT with cPTIO treatment). All experiments were replicated three times. Pollen tubes snap frozen in liquid nitrogen after collection and then stored at −80 °C before use.

RNA extraction, library preparation and RNA-Seq

Nine pollen tube samples from three independent experiments were collected after incubated under control, low temperature and NO treatment. For Illumina sequencing, total RNA was extracted and treated with RNAiso Plus (TaKaRa, Japan) and RNase-free DNase I (Takara Biotechnology, China), respectively, according to the manufacturer’s instructions. We confirmed the integrity, quality and quantity of total RNA using Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and a spectrophotometer (NanoDrop, Wilmington, DE), and used TruSeq RNA Sample Prep Kit v2 (Illumina) to purify 200 ng total RNA for each sample by oligo-dT. Then the total RNA was fragmented using Elute, Prime and fragmentation buffer. We produced the first-strand cDNA using First Strand Master Mix and Super Script II (Invitrogen) and created the second-strand cDNA using Second Strand Master Mix. The double-stranded DNA was then purified by QiaQuik PCR purification kit (Qiagen, Beijing, China) before end repaired by incubation with End Repair Mix at the 30 °C for 30 min. Subsequently, the end-repaired cDNA was added single A base followed by adapter ligation. We performed PCR with PCR Primer Cocktail and PCR Master Mix to amplify the cDNA. The final library was quantified determining the average molecule length measured by the Agilent 2100 bioanalyzer (Agilent DNA 1000 Reagents) and via QPCR (TaqMan Probe), the cDNA library was eventually sequenced on the HiSeq 2000 System (TruSeq SBS KIT-HS V3,Illumina) in the Beijing Genomics Institute (Shenzhen, China; http://www.genomics.cn/index). Above all, the RNA extraction, library preparation and RNA-Seq were carried out according to descriptions showed by Wang et al. [5] and Wang et al. [111].

De novo assembly and data filtering

Raw reads produced directly from the sequencing platform include dirty reads containing unknown or low-quality bases and adapters. After removing the reads with adaptors, the reads with more than 5 % unknown nucleotides and the low-quality reads (>20 % read rate with a quality value ≤ 10), the remaining clean reads were used for the subsequent analysis [112]. Trinity [112] software was then used to assemble short strips. The sequences resulting from Trinity are called unigenes. Using TGICL [113] to further account for splicing and redundancy, these sequences were assembled into homologous transcript clusters. The same process was applied to all three samples to obtain the longest possible non-redundant unigenes. Based on the clustering of the homologous transcripts, the unigenes were divided into clusters composed of several high-similarity (>70 %) unigenes (starting with CL) and singletons (starting with unigenes). Finally, the unigenes were aligned by BLASTX (E-value < 0.00001) to protein databases (in order of priority) NR, Swiss-Prot, KEGG and COG. The proteins showed the highest ranks in the results of BLAST search were used to determine the coding region sequences of the unigenes. We used ESTS software [114] to determine the sequence direction, if a unigene could not align with any of the above databases. Thereby the nucleotide sequence in the direction from 5’ to 3’ and the amino sequence of the predicted coding region are produced. We have deposited all sequencing data at the sequence read archive of NCBI (Accession number SRR3270364, SRR3270376, SRR3270829, SRR3270928, SRR3270974, SRR3270993, SRR3270997, SRR3271001 and SRR3271002).

DEGs selection

We used Tags per Million (TPM) to normalize clean tags of each library to obtain the normalized gene expression levels. DEGs were determined using the NOISeq-bio method as previously reported by Tarazona et al. [115]. We estimated the fold changes (log₂Ratio) on the basis of the normalized gene expression levels in each sample. We used the absolute values of log₂Ratio > 1 and probability > 0.7 as the thresholds criteria for significant differences in gene expression [31]. In addition, DEGs related to plant hormone metabolism and signaling, TFs and the vesicle polarized trafficking, cell wall biosynthesis, the ubiquitination machinery of the ubiquitin system and species-specific secondary metabolites pathways were mainly analyzed.

Unigene annotation and classification

Bioinformatics procedures were employed to classify the annotation of unigenes. BLAST was used against the NT, NR, KEGG, Swiss-Prot and COG (E-value ≤ 1.0E-5) databases for annotation, and subsequently the number of unigenes annotated within each database was counted. Based on the NR annotation, GO functional annotations can be obtained. First, the Blast2GO [32] to get the GO annotation of the unigenes, while the WEGO [33] was used to analyze GO functional classifications for all unigenes and determine the distribution of species gene functions at the macro level. The KEGG annotation was employed to get pathway annotations for our unigenes.

Clustering of gene expression profiles

We performed a hierarchical cluster analysis using the log₂Ratio of the 125 co-expressed DEGs in the CK-VS-LT and CK-VS-NO expression libraries by cluster [116] and Java Treeview [117] software.

qRT-PCR analysis

The qRT-PCR analysis was carried out according to the description shown by Wang et al. [111]. The primers used in this assay were shown in Additional file 16: Table S14, and the C. sinensis 18S rRNA was employed as the reference.

Statistical analyses

Statistical analyses were performed as previous reported by Wang et al. [111]. Briefly, all experiments were replicated at least three times, and all data were represented as the means ± standard deviations (SD). Group differences were tested using one-way ANOVA and Duncan’s test, and significant differences are represented by different letters (P < 0.05). Data analysis was performed using SPSS 20 software.