Background {#Sec1}
==========

Tea (*Camellia sinensis* (L.) O. Kuntze) is a popular evergreen beverage tree, grown worldwide in different agro-climatic between 45°N and 34°S latitude \[[@CR1]\]. Low-temperature stress is one of the most critical abiotic stresses that negatively influence the growth, development and geographical distribution of tea plants \[[@CR2]\]. Previous studies have noted that low temperature lead to a weakened photosynthetic capacity and strengthened respiration rates \[[@CR3]\], and a total of 105 metabolic processes were changed in tea plants, such as carbohydrate metabolism and amino acid metabolism \[[@CR4]\], which ultimately decreased the yield and quality of the tea plants. Thus, finding the mechanisms of improving the tolerance of *C. sinensis* to low temperature is of great importance. Recently, certain molecular mechanisms driving the responses of tea plants to low temperature, such as Ca^2+^ signaling, APETALA2/ethylene responsive factor (AP2/ERF), basic helix-loop-helix (bHLH), WRKY and MYB domain-containing protein (MYB) transcription factors (TFs) and zinc finger protein (ZFP), have been demonstrated, particularly in vegetative tissues \[[@CR5]\]. These findings provided a solid foundation for further exploring the cold-related genes and improved the understanding of plant tolerance to low-temperature stress. However, previous reports revealed that gene expression in reproductive tissues was different from that in vegetative tissues during cold tolerance \[[@CR6]\], implying that reproductive tissue responds to low-temperature stress by different mechanisms \[[@CR7]\]. Furthermore, the pollen tube, as a key component of genetic breeding, plays a crucial role in the higher plants' reproduction, and its growth is highly susceptible to low temperature. For example, certain reports indicate that low temperature clearly reduce pollen tube elongation and disrupt its morphology \[[@CR8], [@CR9]\]. Interestingly, Camellia pollen can germinate *in vitro* even at 3 °C, whereas lily pollen cannot germinate at low temperature \[[@CR10]\], implying higher tolerance to low temperature in Camellia pollen. Our previous report also demonstrated that *C. sinensis* pollen could germinate and grow at low temperature, although the germination and elongation were both retarded by cold stress \[[@CR11]\]. Thus, the basic molecular mechanisms regulating pollen germination and tube elongation in response to low-temperature stress in *C. sinensis* deserve further attention.

Recently, nitric oxide (NO) has gained recognition for its significance in various essential physiological progresses in plants \[[@CR12]\], and much attention has been focused on its roles in the plant tolerance to diverse abiotic stresses, including high salt \[[@CR13]\], drought \[[@CR14]\], heat \[[@CR15]\], heavy metal \[[@CR16]\], ultraviolet radiation \[[@CR17]\], and cold stress \[[@CR18]\]. Evidence suggests that NO production regulates diverse responses during cold stress. For example, NO production was increased in *Triticum aestivum* roots under treatment with cold stress \[[@CR19]\], whereas it was reduced in *Capsicum annuum* leaves under treatment with low temperature \[[@CR20]\]. Furthermore, NO has been confirmed to participate in the anabolic pathway and signal transduction of plant hormones under low-temperature stress. Guo et al. \[[@CR21]\] revealed that NO participates in mediating the cold-induced expression of an enzyme related to the ethylene precursor *MfSAMS1* (S-adenosylmethionine synthetase) by interacting with abscisic acid (ABA) and H~2~O~2~ in *Medicago sativa* subsp. *Falcate*. Recent reports have shown that NADPH oxidase inhibitors and reactive oxygen species scavengers reduced 24-epibrassinolide (EBR)-induced NO production, whereas EBR-induced H~2~O~2~ synthesis was not sensitive to NO. Decreased NO content inhibited EBR-induced low-temperature tolerance and partly inhibited EBR-induced expression of genes coding for several antioxidant enzymes and their relative activities in cucumber, suggesting that NO plays key roles in the H~2~O~2~-dependent induction of plant tolerance to low-temperature stress by brassinosteroids (BRs) \[[@CR22]\]. Additionally, accumulating evidence has suggested that TFs are modulated by NO during cold stress. For instance, the C repeat binding factor 1 (CBF1) and CBF3 TFs are accumulated because of NO during *Arabidopsis* responding to low temperature \[[@CR23]\]. Similarly, application of the NO donor sodium nitroprusside (SNP) increased *LeCBF1* expression, whereas nitro-arginine (a NO synthesis competitive inhibitor) reduced its expression in tomato fruits during cold tolerance \[[@CR24]\]. A recent study has indicated that NO participates in ABA-dependent cold signaling by negatively regulating the DNA binding of MYB2 by S-nitrosylation \[[@CR25]\]. However, NO's participation in low-temperature response is controlled by sophisticated signaling crosstalk for which the potential molecular mechanisms are still undeciphered.

As shown in the previous reports, NO plays crucial roles in pollen tube tip growing, for example its participation in the pollen tube growth regulation and reorientation in lily and *Arabidopsis* \[[@CR26], [@CR27]\]. In addition, extracellular nucleotides suppress *Arabidopsis* pollen germination and tube growth, which is modulated partly through the NO signaling pathway \[[@CR28]\]. Wang et al. \[[@CR29]\] reported that NO affects the cell wall construction in pollen tubes through the alteration of extracellular Ca^2+^ influx and F-actin organization, which consequently affects pollen tube elongation in *Pinus bungeana*. More interestingly, our previous report revealed that NO production from a NO synthase (NOS)-like enzymatic reaction enhanced the inhibiting effect of low-temperature stress on pollen germination and tube elongation in *C. sinensis*, and decreased free proline accumulation, which was accomplished partly through the guanosine 3', 5'-cyclic monophosphate (cGMP) signaling pathway \[[@CR11]\]. However, the downstream signaling and function mechanisms behind the NO-induced reduction of pollen tube elongation under low-temperature stress in *C. sinensis* remain unclear.

High-throughput transcriptomic analysis techniques can provide an overview of expression changes in response to biotic and abiotic stresses. In particular, RNA-Seq has become a widely used method for the analysis of species whose genome sequence is not yet available \[[@CR30]\]. To understand the potential mechanisms driving NO-modulated pollen tube growth during exposure to low-temperature stress, we utilized the Illumina HiSeq™2000 platform to sequence the transcriptomes of *C. sinensis* pollen tubes cultured at 25 °C and low temperature and exposed to exogenous NO treatment, and all of the experiments were replicated three times to guarantee the reliability of the results in this study. Furthermore, the differential gene expression from the transcriptomes of the three treatments was systematically examined, the differentially expressed genes (DEGs) related to plant hormone signaling pathways, TFs, vesicle polarized trafficking, cell wall biosynthesis, the ubiquitination machinery of the ubiquitin system and species-specific secondary metabolite pathways were identified and analyzed. The assembled, annotated transcriptome will present a useful genomic resource for facilitating the investigations on the molecular mechanisms of NO's involvement in pollen tube tolerance to low-temperature stress.

Results {#Sec2}
=======

Illumina sequencing and *de novo* assembly {#Sec3}
------------------------------------------

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^TM^2000 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](#Tab1){ref-type="table"}. 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](#Tab2){ref-type="table"}). These unigenes size distribution was displayed in Additional file [1](#MOESM1){ref-type="media"}: Figure S1. As shown in Fig. [1](#Fig1){ref-type="fig"}, 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.Table 1Summary of the output in the transcriptomes of pollen tubesSampleTotal raw readsTotal clean readsTotal clean nucleotides (nt)Q20 percentageN percentageGC percentageCsPT-CK-154,971,74851,933,9884,674,058,92097.92 %0.00 %45.44 %CsPT-CK-257,928,04654,556,8164,910,113,44097.92 %0.00 %45.23 %CsPT-CK-356,727,45253,686,7564,831,808,04097.93 %0.00 %45.39 %CsPT-LT-158,915,63055,506,6664,995,599,94097.86 %0.00 %45.64 %CsPT-LT-255,468,88052,434,8544,719,136,86097.85 %0.00 %45.39 %CsPT-LT-354,793,79451,330,9084,619,781,72097.77 %0.01 %45.35 %CsPT-NO-153,553,77250,257,9764,523,217,84097.84 %0.01 %45.20 %CsPT-NO-255,260,05851,851,1064,666,599,54097.82 %0.01 %45.40 %CsPT-NO-354,754,89851,364,4664,622,801,94097.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 Table 2Summary of the assembly quality of the transcriptomes of the pollen tubesSampleTotal numberTotal length(nt)Mean length(nt)N50Total consensus sequencesDistinct clustersDistinct singletonsContigCsPT-CK-160,61120,192,162333571CsPT-CK-261,24620,491,956335581CsPT-CK-363,11720,892,307331561CsPT-LT-156,10219,226,389343599CsPT-LT-258,62219,960,405340595CsPT-LT-360,54720,132,613333572CsPT-NO-159,75019,797,938331559CsPT-NO-258,88419,734,167335575CsPT-NO-359,55019,805,943333568UnigeneCsPT-CK-139,72625,741,808648116539,72614,09525,631CsPT-CK-240,44026,270,780650115240,44014,71725,723CsPT-CK-341,62627,483,327660120441,62615,31026,316CsPT-LT-136,99324,057,930650118236,99313,29023,703CsPT-LT-239,07025,937,230664120339,07014,22424,846CsPT-LT-339,43925,906,326657119139,43914,30125,138CsPT-NO-139,51425,545,806647117339,51414,21125,303CsPT-NO-238,29824,889,732650116938,29813,67624,622CsPT-NO-339,06124,971,317639115339,06113,92325,138All45,43246,435,0791022157445,43222,81222,620Total 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 Fig. 1Pearson correlation analysis of the replicates from the control (CK), low-temperature and NO treatments on the *Camellia sinensis* pollen tubes. **a**, **b** and **c** are analyses of the CK treatment; **d**, **e** and **f** are analyses of the low-temperature (LT) treatment, and **g**, **h** and **i** are analyses of the NO donor DEA NONOate (NO) treatment.

Annotation of predicted proteins {#Sec4}
--------------------------------

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](#Fig2){ref-type="fig"}, 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^−60^), 20.4 % showed a strong homology (10^−60^ \< E-value \< 10^−30^), and an additional 31.3 % showed homology (10^−30^ \< E-value \< 10^−5^) (Fig. [2a](#Fig2){ref-type="fig"}). 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](#Fig2){ref-type="fig"}). 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 balsamifera*s ubsp*. trichocarpa* (6.7 %) (Fig. [2c](#Fig2){ref-type="fig"}).Fig. 2Homology search of unigenes against the NR database. **a** and **b** show the E-value and similarity distribution of the top BLAST hits for each unique sequence, respectively; **c** shows species distribution of the top BLAST hits for every homologous sequence

Functional classification by COG {#Sec5}
--------------------------------

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](#Fig3){ref-type="fig"}). 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".Fig. 3COG functional classification of the transcriptomes. Among 33,716 hits in the NR database, 13,029 unigenes with significant homologies in the COG database (E-value \< 0.00001) were classified into 24 COG categories

DEGs and GO classification {#Sec6}
--------------------------

We filtered the DEGs (Additional file [2](#MOESM2){ref-type="media"}: Table S1, Additional file [3](#MOESM3){ref-type="media"}: Table S2 and Additional file [4](#MOESM4){ref-type="media"}: Table S3) using the absolute values of log~2~Ratio \> 1 and probability \> 0.7 as the thresholds according to Zhao et al. \[[@CR31]\]. 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](#Fig4){ref-type="fig"}). 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](#Fig4){ref-type="fig"}).Fig. 4The number of DEGs identified in the comparisons of CK-VS-LT, CK-VS-NO and LT-VS-NO. CK: control; LT: low-temperature treatment; NO: NO donor DEA NONOate treatment. CK-VS-LT: comparison between CK and LT. CK-VS-NO: comparison between CK and NO. LT-VS-NO: comparison between LT and NO. The absolute values of log~2~Ratio \> 1 and probability \> 0.7 were used as threshold for assigning significance. **a** the number of up-regulated (*blue*) and down-regulated (*red*) DEGs in the comparisons of CK-VS-LT, CK-VS-NO and LT-VS-NO are shown, respectively. **b** the number of specifically and co-expressed DEGs from the comparisons of CK-VS-LT, CK-VS-NO and LT-VS-NO are shown, respectively

GO classifies genes on the basis of three functional ontologies, that is, Molecular function, Cellular component and Biological process. We used Blast2GO \[[@CR32]\] to obtain GO annotations and Web Gene Ontology Annotation Plot (WEGO) \[[@CR33]\] 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](#MOESM5){ref-type="media"}: Table S4, Additional file [6](#MOESM6){ref-type="media"}: Table S5 and Additional file [7](#MOESM7){ref-type="media"}: 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](#Fig5){ref-type="fig"}). As shown in Fig. [5](#Fig5){ref-type="fig"}, 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.Fig. 5GO functional classifications of DEGs found from the comparisons of CK-VS-LT, CK-VS-NO and LT-VS-NO are annotated in the categories of Biological process (*blue*), Cellular component (*red*) and Molecular function (*green*). **a** GO classifications of DEGs from comparison of CK-VS-LT; **b**: GO classifications of DEGs from comparison of CK-VS-NO; **c**: GO classifications of DEGs from comparison of LT-VS-NO

Cluster analysis of DEGs co-expressed in CK-VS-LT and CK-VS-NO {#Sec7}
--------------------------------------------------------------

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](#MOESM8){ref-type="media"}: 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](#Fig6){ref-type="fig"}). The largest group (Group 1) contained 71 genes (56.8 %) (from Unigene 11102_All to Unigene 20419_All) (Additional file [8](#MOESM8){ref-type="media"}: 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.Fig. 6Clustering analysis of the DEGs co-expressed in CK-VS-LT and CK-VS-NO. This analysis is based on the log~2~Ratio values of DEGs. Columns displays the two comparisons (CK-VS-LT and CK-VS-NO), and rows represent the 125 DEGs. Up- and down-regulated gene expressions are showed by *red* and *green*, respectively. The 125 genes exhibit three expression patterns

DEGs identified from pollen tubes treated with low temperature and NO {#Sec8}
---------------------------------------------------------------------

In this study, numerous genes were observed that exhibited changed expression levels when pollen tubes were exposed to low temperature or NO. These genes included factors related to plant hormone metabolism and signaling pathways, such as ABA8ox and type 2C protein phosphatase (PP2C) which participate in the ABA metabolism and signaling network (Tables [3](#Tab3){ref-type="table"} and [4](#Tab4){ref-type="table"}); TFs, such as MYB and ZFP transcription factors (Tables [3](#Tab3){ref-type="table"} and [4](#Tab4){ref-type="table"}); vesicle polarized trafficking and cell wall biosynthesis (Additional file [9](#MOESM9){ref-type="media"}: Table S8 and Additional file [10](#MOESM10){ref-type="media"}: Table S9); the ubiquitination machinery of the ubiquitin system (Additional file [11](#MOESM11){ref-type="media"}: Table S10 and Additional file [12](#MOESM12){ref-type="media"}: Table S11); and species-specific secondary metabolite pathways (Additional file [13](#MOESM13){ref-type="media"}: Table S12 and Additional file [14](#MOESM14){ref-type="media"}: Table S13). Furthermore, many of these DEGs responded to both low temperature and NO.Table 3DEGs related to plant hormone metabolism and signaling pathways and transcription factors (TFs) between CK and LT (CK-VS-LT)GeneIDGene lengthlog~2~Ratio(LT/CK)Up-Down-Regulation(LT/CK)ProbabilityGene annotationUnigene14044_All256−1.843769497down0.716541745abscisic acid 8\'-hydroxylaseCL580.Contig48_All7541.27259979up0.774663281PP2CUnigene9712_All720−1.700722725down0.762634179brassinosteroid insensitive 1Unigene4038_All663−1.15134991down0.740396494brassinosteroid insensitive 1CL4145.Contig1_All436−1.824955819down0.800581069BRs 1-associated receptor kinaseUnigene4486_All236−2.751491338down0.84984985Ethylene response sensor 1Unigene22534_All230−1.419906337down0.738457156Ethylene response sensor 1CL2610.Contig1_All1080−1.518999731down0.801879115EIN3-binding F-box proteinUnigene10222_All255−2.230648476down0.7352979EIN3-binding F-box proteinCL2941.Contig1_All5431.596587443up0.813476078auxin-repressed proteinCL2941.Contig2_All10621.516730948up0.810269619auxin-repressed proteinUnigene7793_All2381.145845871up0.773991064Auxin response factorUnigene13496_All3791.090218234up0.761820357indole-3-acetic acid-amido synthetaseUnigene21651_All201−1.631827096down0.734695671indole-3-acetic acid-induced proteinUnigene15959_All391−1.60647149down0.739220521DELLA proteinUnigene10398_All4561.173251983up0.732937816DELLA proteinUnigene11055_All1276−1.283576382down0.760756692Gibberellin 2 oxidaseUnigene1092_All683−1.991364704down0.737188408Gibberellin 20 oxidaseUnigene9129_All9891.748223423up0.822590884Gibberellin receptor GID1Unigene21715_All207−5.94952025down0.749639884MYB transcription factorUnigene17702_All299−1.640960381down0.738882785MYB transcription factorUnigene5788_All609−2.816635393down0.906141914Zinc finger proteinUnigene21920_All236−2.145887308down0.748960342Zinc finger proteinUnigene9475_All3821.147076601up0.709884681Zinc finger proteinCL1123.Contig1_All1049−1.093093550.716326082downMADS-box transcription factorUnigene19750_All3641.6834475470.796843185upMADS-box transcription factorUnigene721_All691−1.0999646250.762959708downMADS-box transcription factorUnigene7926_All4001.0649328470.711558713upMADS-box transcription factorCL1553.Contig1_All531−1.0269238680.758009229downMADS-box transcription factorUnigene2215_All492−1.992582568down0.704653434CAMTAUnigene1182_All690−4.091560517down0.945819804ERF transcription factorUnigene8478_All2571.323055163up0.787730006ERF transcription factorUnigene7269_All1007−1.603627707down0.77292333WRKY transcription factorUnigene22002_All292−2.6131679760.704958617downbHLH transcription factorUnigene1490_All285−1.7742052670.79511951downbHLH transcription factorThe absolute values of the log~2~Ratio (LT/CK) \> 1 and probability \> 0.7 were used as threshold for assigning significance. *CK* control, *LT* 4 °C treatment Table 4DEGs involved in the plant hormone metabolism and signaling pathways and transcription factors (TFs) between CK and NO (CK-VS-NO)GeneIDGene lengthlog~2~Ratio(NO/CK)Up-Down-Regulation(NO/CK)ProbabilityGene annotationUnigene14044_All256−1.843769497down0.716541745abscisic acid 8\'-hydroxylaseCL580.Contig48_All7541.070746844up0.727803027PP2CUnigene4038_All663−1.137884213down0.84974181Leucine-rich repeat (LRR) proteinUnigene22534_All230−1.154935001down0.706870123Ethylene response sensor 1Unigene17743_All2721.415347606up0.7430468indole-3-acetic acid-amido synthetaseUnigene21651_All201−1.525739727down0.723650002indole-3-acetic acid-induced proteinUnigene17702_All299−1.618842865down0.747241362MYB transcription factorUnigene4344_All557−1.108187448down0.725592533MYB transcription factorUnigene21920_All236−1.849780074down0.72892945Zinc finger proteinCL1553.Contig1_All531−1.28555up0.738229MADS-box transcription factorUnigene1516_All2301.242887153up0.792390331CAMTAUnigene2215_All492−1.171631887down0.701307237CAMTAUnigene20631_All228−5.122748755down0.78094284ERF transcription factorUnigene13943_All358−1.523794549down0.785566632bHLH transcription factorCL898.Contig1_All231−2.128187976down0.744496571bHLH transcription factorThe absolute values of log~2~Ratio (NO/CK) \> 1 and probability \> 0.7 were used as threshold for assigning significance. *CK* control, *NO* NO treatment

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

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^2+^ signaling and TFs played an important role in this process \[[@CR11], [@CR25]\]. 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 \[[@CR34]\] and PLD activation in *C. bungeana* exposed to low-temperature stress was related to a drop in Ca^2+^ content in membrane fractions \[[@CR35]\]. Moreover, TOPLESS (TPL) was involved in mutual adjustment between hormones via NO in *Petunia hybrida* \[[@CR36]\]. Additionally, genes encoding a non-specific serine/threonine protein kinase were down-regulated expressed during response to low temperature in two grapevine varieties \[[@CR37]\]. 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 \[[@CR38]\]. 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](#Fig7){ref-type="fig"}). 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.Fig. 7QRT-PCR analysis confirms the expression patterns indicated by the sequencing. Nine genes were randomly selected for qRT-PCR analysis and their primers were presented in Additional file [16](#MOESM16){ref-type="media"}: Table S14. The *C. sinensis* 18S rRNA was employed as reference

Discussion {#Sec10}
==========

According to our previous report, NO negatively regulated *C. sinensis* pollen tube growth under low-temperature stress partially through the cGMP signaling pathway \[[@CR11]\]. However, the potential molecular mechanism during this process is still unknown. Here, we comprehensively analyzed the gene expression changes in *C. sinensis* pollen tubes treated with low temperature and exogenous NO by RNA-Seq technology to investigate the pollen tube's defense response to low temperature. Genes related to metabolism and signaling pathways of plant hormones, TFs, vesicle polarized trafficking and cell wall biosynthesis, the ubiquitination machinery of the ubiquitin system and species-specific secondary metabolite pathways were the main focus of this study.

Comparison of unigenes between *C. sinensis* leaves and pollen tubes {#Sec11}
--------------------------------------------------------------------

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 \[[@CR5]\]. 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 \[[@CR5]\]. Among 71,289 transcripts based on NR annotation, 30,115 transcripts matched 80,176 GO term annotations \[[@CR5]\]. 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](#MOESM15){ref-type="media"}: 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](#MOESM15){ref-type="media"}: 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](#MOESM15){ref-type="media"}: Figure S2B). These differences between transcriptomes from leaves and pollen tubes might result from the diversity of reproductive and vegetative tissues of plant \[[@CR6]\], further confirming that reproductive tissue responds to low temperature stress using distinct mechanisms from vegetative tissues \[[@CR7]\].

Plant hormone metabolism and signaling pathways {#Sec12}
-----------------------------------------------

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 \[[@CR39]\], and previous investigations have demonstrated that NO enhances tolerance to drought stress via promoting ABA-induced stomata closure \[[@CR40], [@CR41]\]. In our study, we identified one *ABA8ox* gene (Unigene14044_All) and one *PP2C* gene (CL580.Contig48_All) in the CK-VS-LT (Table [3](#Tab3){ref-type="table"}) and CK-VS-NO libraries (Table [4](#Tab4){ref-type="table"}). As shown in previous reports, NO-induced ABA content reduction in *Arabidopsis* correlates with the modulation of ABA8ox protein expression \[[@CR42]\]. 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 \[[@CR43]\]. As a negative regulator, PP2C inactivates SnRK2, thereby regulating ABA signals \[[@CR44]\]. 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 \[[@CR45]\]. 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 \[[@CR46]\]. 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. \[[@CR47]\]. We identified two *BRI1* genes (Unigene9712_All and Unigene4038_All) and one *BAK1* gene (CL4145.Contig1_All) in CK-VS-LT (Table [3](#Tab3){ref-type="table"}) as well as one *BRI1* gene (Unigene4038_All) in CK-VS-NO (Table [4](#Tab4){ref-type="table"}). 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* \[[@CR48]\]. 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 \[[@CR49]\]. BR-induced NO production and NO-activated ABA biosynthesis are significant mechanisms for BR-increased tolerance to water-stress \[[@CR50]\]. 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 \[[@CR51], [@CR52]\]. Previous reports indicated that ET signaling negatively regulates *Arabidopsis* tolerance to cold stress by inhibiting the expression of *CBF* and *type-A ARR* genes \[[@CR53]\]. Ethylene response sensor 1 (ERS1) has been identified as one of the membrane-located receptor proteins in the ET signaling pathway \[[@CR54]\]. 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* \[[@CR55]\]. 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 \[[@CR56]\]. EIN3-Binding F-box-1 (EBF1) and EBF2 tightly regulate the stability of EIN3 at the protein level \[[@CR57]\], and EBF2 transcription is in turn activated by EIN3, thereby forming a negative feedback loop \[[@CR58]\]. Interestingly, recent work has shown that ERS1 and EBF1/2 are differentially expressed in Papaya fruit ripening disorders caused by cold injury \[[@CR59]\]. 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](#Tab3){ref-type="table"}). 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 \[[@CR60]\], salt stress \[[@CR61]\], UV-B stress \[[@CR62]\] and hypersensitive responses \[[@CR63]\]. One *ERS1* gene (Unigene22534_All) was also identified in the CK-VS-NO library (Table [4](#Tab4){ref-type="table"}). 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 \[[@CR64], [@CR65]\]. 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 \[[@CR66]\]. 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](#Tab3){ref-type="table"}). Interestingly, as shown in a recent report, auxin participates in regulating pollen tube tip growth \[[@CR67]\]. Therefore, it is reasonable to speculate that auxin intervenes in pollen tube responding to low-temperature stress in *C. sinensis*. Furthermore, Markus et al. \[[@CR68]\] 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](#Tab4){ref-type="table"}). 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 \[[@CR69]\]. 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 \[[@CR70], [@CR71]\]. 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 \[[@CR72]\]. 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](#Tab3){ref-type="table"}). 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. \[[@CR73]\] 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) {#Sec13}
---------------------------

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 \[[@CR74]\]. Similarly, the soybean typical Cys2/His2-type (C2H2-type) zinc finger *GmZF1* enhances cold stress tolerance by regulating gene expression in transgenic *Arabidopsis* \[[@CR75]\]. In our study, two *MYB* genes and three *ZFP* genes were identified in CK-VS-LT as shown in Table [3](#Tab3){ref-type="table"}. 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 \[[@CR76]\]. 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 \[[@CR77]\]. 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](#Tab3){ref-type="table"} and [4](#Tab4){ref-type="table"}), 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 \[[@CR25]\]. 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* \[[@CR78]\]. 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 \[[@CR79]--[@CR81]\]. The CaM-interacting protein group CAMTA has been identified in the plant defense signaling cascade \[[@CR82]\]. Moreover, MADS-box \[[@CR83]\] and CAMTA \[[@CR84]\] 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](#Tab3){ref-type="table"}). 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](#Tab3){ref-type="table"} and [4](#Tab4){ref-type="table"}). Previous experiments have demonstrated that NO modulates Ca^2+^ signaling in lily and *Arabidopsis* pollen tubes through directly imaging the concentration of cytosolic free Ca^2+^ during NO-induced tube re-orientation \[[@CR26], [@CR27]\]. We speculate that the *CAMTA* gene (Unigene2215_All) might participate in the regulation of Ca^2+^ 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 \[[@CR85]\], and the metabolism of arginine is regulated by MADS-box \[[@CR86]\]. 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](#Tab3){ref-type="table"}) as well as one *ERF* gene and two *bHLH* genes in CK-VS-NO (Table [4](#Tab4){ref-type="table"}). Using forward and reverse genetics, the ERF \[[@CR51]\], WRKY \[[@CR87]\] and bHLH \[[@CR88]\] 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 {#Sec14}
----------------------------------------------------------------------

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 \[[@CR89]\]. 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 \[[@CR90]\]. Furthermore, PI 4-phosphate 5-kinase (PIP5K) isoforms can catalyze the PI 4, 5-bisphosphate (PIP2) synthesis, which eventually regulates vesicle trafficking \[[@CR91]\]. Additionally, Feng et al. \[[@CR92]\] 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* \[[@CR93]\]. Here, we identified four *PI4K* genes, one *PIP5K* gene and two *VAMP* genes in CK-VS-LT (Additional file [9](#MOESM9){ref-type="media"}: 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](#MOESM10){ref-type="media"}: Table S9). Interestingly, NO involvement in pollen tube growth regulation has been reported to occur via the modulation of vesicle polarized trafficking \[[@CR29]\], and NO is also required for vesicle formation and trafficking in root hair tip growth in *Arabidopsis* \[[@CR94]\]. 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 \[[@CR95]\] and could be regulated by NO in *Pinus bungeana* \[[@CR29]\]. COBRA-like genes are important factors for secondary cell wall development in the root growth zone of *Arabidopsis* \[[@CR96]\]. We also identified one *COBRA-like* gene (Unigene10170_All) in the DEGs between CK-VS-LT and CK-VS-NO (Additional file [9](#MOESM9){ref-type="media"}: Table S8 and Additional file [10](#MOESM10){ref-type="media"}: 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 {#Sec15}
------------------------------------------------

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 \[[@CR97]\]. 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 \[[@CR98]\] 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 \[[@CR99]\]. Additionally, CRLs can consist of a substrate-recruiting protein, such as F-box \[[@CR100]\] and broad complex, tramtrack, bric-a-brac/pox virus and zinc finger (BTB/POZ) \[[@CR101]\]. 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](#MOESM11){ref-type="media"}: 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 \[[@CR102]\], 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 \[[@CR98]\]. Previous reports have shown that multiple E3s regulate ABA-dependent stress signaling \[[@CR103]\]. 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* \[[@CR104]\], 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](#MOESM12){ref-type="media"}: 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 {#Sec16}
-----------------------------------------------

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 \[[@CR105]\]. Furthermore, the accumulation of catechins in *C. sinensis* is very susceptible to low temperature \[[@CR106]\]. 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 \[[@CR107]\]. 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 \[[@CR108]\]. 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* \[[@CR109], [@CR110]\]. In the present study, some genes related to flavonoids, caffeine and theanine biosynthesis pathways were down-regulated in CK-VS-LT (Additional file [13](#MOESM13){ref-type="media"}: Table S12). In addition, we also identified some genes in CK-VS-NO (Additional file [14](#MOESM14){ref-type="media"}: 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).

Conclusions {#Sec17}
===========

We produced a dataset containing 45,432 unigenes from the *C. sinensis* pollen tube transcriptomes using paired-end and *de novo* sequencing with the Illumina HiSeq™2000 platform. Among these unigenes, 36,097 were annotated with descriptions from the NR, NT, Swiss-Prot, KEGG, COG and GO databases. Our dataset includes comprehensive analysis in sequence and DEGs profiling data that provide a dynamic perspective on transcriptomic variations caused by NO and low-temperature in *C. sinensis* pollen tube elongation. Many DEGs were simultaneously found in both the CK-VS-LT and CK-VS-NO libraries in different functional pathways, such as plant hormone signaling (ABA, BRs, ET and auxin), TFs (MYB, ZFP, CAMTA and MADS-box), vesicle polarized trafficking and cell wall biosynthesis (PI4K and COBRA-like), the ubiquitination machinery of the ubiquitin system (E3s, 26S proteasome, U-box and BTB/POZ) and caffeine and theanine biosynthesis (caffeine synthase and glutamate receptor). Combining our findings and previous reports, we propose a crosstalk network which may indicate the probable mechanisms underlying NO mediating *C. sinensis* pollen tube growth during the response to low-temperature stress (Fig. [8](#Fig8){ref-type="fig"}) \[[@CR102], [@CR111]\]. Plant hormone signaling pathways (ABA, BRs, ET and auxin) respond to pollen germination and tube elongation at low temperature, and then TFs (MYB, ZFP, CAMTA and MADS-box) modulated relevant proteins expression. Moreover, the UPS (E3s, 26S proteasome, U-box and BTB/POZ) probably involves in affecting pollen germination and tube growth in response to low-temperature stress by modulating the activity, localization or stability of different proteins. Additionally, Ca^2+^ gradient, vesicle polarized trafficking as well as cell wall biosynthesis might participate in pollen tube elongation under low-temperature stress through the NO signaling pathway by CAMTA, PI4K and COBRA-like. Taken together, these findings offer a more profound understanding of the candidate genes employed by *C. sinensis* pollen tubes to integrate low temperature signals and finely tune their tolerance responses, present a fundamental transcriptomic resource to further explore the molecular mechanisms of reproductive tissue tolerance to low-temperature.Fig. 8Model of the potential signaling pathway of NO in *C. sinensis* pollen tube response to low temperature. This summarized model was based on the results reported by Wang et al. \[[@CR111]\] and Lyzenga *ea al*. \[[@CR102]\]. Plant hormone signaling pathways (ABA, BRs, ET and auxin) respond to pollen germination and tube elongation at low temperature, and then TFs (MYB, ZFP, CAMTA and MADS-box) modulated relevant proteins expression. Moreover, the UPS (E3s, 26S proteasome, U-box and BTB/POZ) probably involves in affecting pollen germination and tube growth in response to low-temperature stress by modulating the activity, localization or stability of different proteins. Additionally, Ca^2+^ gradient, vesicle polarized trafficking as well as cell wall biosynthesis might participate in pollen tube elongation under low-temperature stress through the NO signaling pathway by CAMTA, PI4K and COBRA-like

Methods {#Sec18}
=======

Plant materials {#Sec19}
---------------

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~3~)~2~  **·** 4H~2~O, 0.01 % H~3~BO~3~, 5 % PEG 4000 and 30 mM MES, pH 6.0) at 25 °C in darkness for 30 min \[[@CR111]\]. 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 {#Sec20}
-----------------------------------------------

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. \[[@CR5]\] and Wang et al. \[[@CR111]\].

*De novo* assembly and data filtering {#Sec21}
-------------------------------------

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 \[[@CR112]\]. Trinity \[[@CR112]\] software was then used to assemble short strips. The sequences resulting from Trinity are called unigenes. Using TGICL \[[@CR113]\] 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 \[[@CR114]\] 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 {#Sec22}
--------------

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. \[[@CR115]\]. We estimated the fold changes (log~2~Ratio) on the basis of the normalized gene expression levels in each sample. We used the absolute values of log~2~Ratio \> 1 and probability \> 0.7 as the thresholds criteria for significant differences in gene expression \[[@CR31]\]. 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 {#Sec23}
-------------------------------------

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 \[[@CR32]\] to get the GO annotation of the unigenes, while the WEGO \[[@CR33]\] 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 {#Sec24}
--------------------------------------

We performed a hierarchical cluster analysis using the log~2~Ratio of the 125 co-expressed DEGs in the CK-VS-LT and CK-VS-NO expression libraries by cluster \[[@CR116]\] and Java Treeview \[[@CR117]\] software.

qRT-PCR analysis {#Sec25}
----------------

The qRT-PCR analysis was carried out according to the description shown by Wang et al. \[[@CR111]\]. The primers used in this assay were shown in Additional file [16](#MOESM16){ref-type="media"}: Table S14, and the *C. sinensis* 18S rRNA was employed as the reference.

Statistical analyses {#Sec26}
--------------------

Statistical analyses were performed as previous reported by Wang et al. \[[@CR111]\]. 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.

Additional files {#Sec27}
================

Additional file 1: Figure S1.Assessment of assembly quality. Distribution of mapped reads within the assembled unigenes determined unigene assembly quality. (TIF 2384 kb) Additional file 2: Table S1.DEGs identified from the comparison between control (CsPT-CK) and 4 °C-treated (CsPT-LT) pollen tbues. All of the samples were replicated three times. CK and LT FPKM: fragments per kb per million reads for each unigene in the CK and LT libraries, respectively. The log~2~Ratio (LT/CK): ratio between the FPKM of LT and CK. The absolute values of log~2~Ratio \> 1 and probability \> 0.7 were used as threshold for assigning significance. Annotation of DEGs against NR, NT, Swiss-Prot protein, KEGG, COG and GO were all reported in the tables. "-": no hit. (XLS 381 kb) Additional file 3: Table S2.DEGs identified from the comparison between control (CsPT-CK) and NO-treated (CsPT-NO) pollen tubes. All of the samples were replicated three times. CK and NO FPKM: fragments per kb per million reads for each unigene in the CK and NO libraries, respectively. The log~2~Ratio (NO/CK): ratio between the FPKM of NO and CK. The absolute values of log~2~Ratio \> 1 and probability \> 0.7 were used as threshold for assigning significance. Annotation of DEGs against NR, NT, Swiss-Prot protein, KEGG, COG and GO were all reported in the tables. "-": no hit. (XLS 250 kb) Additional file 4: Table S3.DEGs identified from the comparison between 4 °C- (CsPT-LT) and NO-treated (CsPT-NO) pollen tubes. All of the samples were replicated three times. LT and NO FPKM: fragments per kb per million reads for each unigene in the LT and NO libraries, respectively. The log~2~Ratio (LT/NO): ratio between the FPKM of LT and CK. The absolute values of log~2~Ratio \> 1 and probability \> 0.7 were used as threshold for assigning significance. Annotation of DEGs against NR, NT, Swiss-Prot protein, KEGG, COG and GO were all reported in the tables. "-": no hit. (XLS 451 kb) Additional file 5: Table S4.GO classification of DEGs identified from the comparison between control (CK) and 4 °C-treated (LT) pollen tubes. (XLS 52 kb) Additional file 6: Table S5.GO classification of DEGs identified from the comparison between control (CK) and NO-treated (NO) pollen tubes. (XLS 40 kb) Additional file 7: Table S6.GO classification of DEGs identified from the comparison between pollen tubes exposed to 4 °C (LT) and pollen tubes treated with NO donor (NO). (XLS 59 kb) Additional file 8: Table S7.The 125 DEGs identified from the comparison between CK-VS-LT and CK-VS-NO. (XLS 89 kb) Additional file 9: Table S8.DEGs involved in vesicle polarized trafficking and cell wall biosynthesis between CK and LT (CK-VS-LT). The absolute values of log~2~Ratio (LT/CK) \> 1 and probability \> 0.7 were used as threshold for assigning significance. CK: control; LT: 4 °C treatment. (DOC 37 kb) Additional file 10: Table S9.DEGs involved in vesicle polarized trafficking and cell wall biosynthesis between CK and NO (CK-VS-NO). The absolute values of log~2~Ratio (NO/CK) \> 1 and probability \> 0.7 were used as threshold for assigning significance. CK: control; NO: NO treatment. (DOC 31 kb) Additional file 11: Table S10.DEGs involved in the ubiquitination machinery of the ubiquitin system between CK and LT (CK-VS-LT). The absolute values of log~2~Ratio (LT/CK) \> 1 and probability \> 0.7 were used as threshold for assigning significance. CK: control; LT: 4 °C treatment. (DOC 55 kb) Additional file 12: Table S11.DEGs involved in the ubiquitination machinery of the ubiquitin system between CK and NO (CK-VS-NO). The absolute values of log~2~Ratio (NO/CK) \> 1 and probability \> 0.7 were used as threshold for assigning significance. CK: control; NO: NO treatment. (DOC 37 kb) Additional file 13: Table S12.DEGs involved in flavonoids, caffeine and theanine biosynthesis pathways between CK and LT (CK-VS-LT). The absolute values of log~2~Ratio (LT/CK) \> 1 and probability \> 0.7 were used as threshold for assigning significance. CK: control; LT: 4 °C treatment. (DOC 36 kb) Additional file 14: Table S13.DEGs involved in flavonoids, caffeine and theanine biosynthesis pathways between CK and NO (CK-VS-NO). The absolute values of log~2~Ratio (NO/CK) \> 1 and probability \> 0.7 were used as threshold for assigning significance. CK: control; NO: NO treatment. (DOC 34 kb) Additional file 15: Figure S2.GO functional classifications in transcriptome data of *C. sinensis* leaves and pollen tubes under low temperature stress. A: GO functional classifications in transcriptome data of *C. sinensis* leaves under cold stress by Wang et al. \[[@CR5]\]. B: GO functional classifications in transcriptome data of *C. sinensis* pollen tubes under cold stress. (TIF 3893 kb) Additional file 16: Table S14.Primers used in the quantitative real-time PCR (qRT-PCR) analysis of DEGs. (DOC 33 kb)

ABA

:   Abscisic acid

ABA8ox

:   ABA 8'-hydroxylase

AP2/ERF

:   APETALA2/ethylene responsive factor

BAK1

:   BRI1-associated kinase 1

bHLH

:   Basic helix-loop-helix

BRI1

:   Brassinosteroid-insensitive 1

BRs

:   Brassinosteroids

BTB/POZ

:   Broad complex tramtrack, bric-a-brac/pox virus and zinc finger

CaM

:   Calmodulin

CAMTA

:   CaM-binding transcription activator

CBF

:   C repeat binding factor

cGMP

:   Guanosine 3' 5'-cyclic monophosphate

COG

:   Cluster of Orthologous Groups of proteins database

cPTIO

:   2-(4-carboxyphenyl)-4,4,5,5-tetramethylimidazoline-1-oxyl-3-oxide

CRLs

:   Cullin--RING ligases

CTR1

:   Constitutive Triple Response 1

DEGs

:   Differentially expressed genes

E1s

:   Ubiquitin activating enzymes

E2s

:   Ubiquitin conjugating enzymes

E3s

:   Ubiquitin ligases

EBF

:   EIN3-Binding F-box

EBR

:   Epibrassinolide

EIL1

:   EIN3-like1

EIN

:   Ethylene insensitive

ERS1

:   Ethylene response sensor 1

ET

:   Ethylene

GA20x

:   GA 20oxidase

GA2ox

:   GA 2-oxidase

GA3x

:   GA 3oxidase

GAs

:   Gibberellins

GID1

:   GA insensitive dwarf 1

GO

:   Gene Ontology database

HECT

:   Homology to E6-Associated Carboxyl-Terminus

IAA

:   Indole-3-acetic acid

KEGG

:   Kyoto Encyclopedia of Genes and Genomesdatabase

LRR

:   Leucine-rich repeat

MYB

:   MYB domain-containing protein

NO

:   Nitric oxide

NOS

:   NO synthase

NR

:   NCBI non-redundant database

NT

:   Mysql-nt database

PHYA

:   Phytochrome A

PHYB

:   Phytochrome B

PI4K

:   Phosphatidylinositol-4-kinase

PI4P

:   Phosphatidylinositol-4-phosphate

PIN1

:   PIN-FORMED 1

PIP2

:   PI 4,5-bisphosphate

PIP5K

:   PI 4-phosphate 5-kinase

PLD

:   Phospholipase D

PP2C

:   Type 2C protein phosphatase

qRT-PCR

:   Real-time quantitative polymerase chain reaction

RING

:   Really Interesting New Gene

RLK

:   LRR receptor-like kinase

SAMS

:   S-adenosylmethionine synthetase

SNF1

:   S-nitrosylation of sucrose non-fermenting 1

SNP

:   Sodium nitroprusside

SnRK

:   SNF1-related protein kinase

TAIR

:   The Arabidopsis Information Resource

TFs

:   Transcription factors

TPL

:   TOPLESS

TPM

:   Tags per Million

UPS

:   Ubiquitin--26S proteasome system

VAMP726

:   Vesicle-associated membrane protein 726

WEGO

:   Web Gene Ontology Annotation Plot

ZFP

:   Zinc finger protein

Not applicable.

Funding {#FPar1}
=======

This work was supported by the National Natural Science Foundation of China (No. 31370014), the Opening Project of State Key Laboratory of Tea Plant Biology and Utilization (No. SKLTOF20150113) and the Fundamental Research Funds For the Central Universities (No. KYZ201509).

Availability of data and materials {#FPar2}
==================================

RNA-Seq data used in this study have been deposited into the NCBI Sequence Read Archive (SRA, <http://www.ncbi.nlm.nih.gov/sra/>) under accession number SRR3270364, SRR3270376, SRR3270829, SRR3270928, SRR3270974, SRR3270993, SRR3270997, SRR3271001 and SRR3271002.

Authors' contributions {#FPar3}
======================

JP, YW and WW conceived the study. JP, WW carried out bioinformatics data analysis. DL, ZS and XY performed qPCR analysis and participated in the sequence alignment. JP and YW wrote the manuscript. YW, WW and PC critically reviewed the manuscript. All authors read and approved the final manuscript.

Competing interest {#FPar4}
==================

The authors declare that they have no competing interests.

Consent for publication {#FPar5}
=======================

Not applicable.

Ethics approval and consent to participate {#FPar6}
==========================================

Not applicable.
