Watermelon (C. lanatus L., cv. Y134) seeds were surface sterilized with 5% sodium hypochlorite (NaOCl) solution for 5 min and then rinsed with running water and distilled water. The sterilized seeds were soaked in distilled water for 8 h and sown directly in pots filled with a mixture of peat/vermiculite (3/1, v/v). Plants were grown in growth chambers with the following environmental conditions: a constant relative humidity of 60–70%, a 12-h photoperiod, 25/18°C (day/night), and a photosynthetic photon flux density of 600 μmol m^-2 s^-1. The plants were watered daily and fertilized with Hoagland’s nutrition solution at 1-day interval.

Seedlings at the four-leaf stage were sprayed with 150 μM melatonin solution for 3 days, with distilled water used as the control. The melatonin (Sigma-Aldrich, St. Louis, MO, USA) solutions were prepared by dissolving the solute in ethanol followed by dilution with Milli-Q water [ethanol/water (v/v) = 1/10,000]. Each plant was sprayed with 20 mL of solution. Twelve hours after the third spray of melatonin, the plants were exposed to cold stress (i.e., 4°C) for 36 h with a 12-h photoperiod and photosynthetic photon flux density of 600 μmol m^-2 s^-1. Leaf samples were harvested at different time-points after imposition of cold stress, such as at 0, 3, 6, 12, 24, and 36 h to analyze Cla020078 (CBF1) and Cla020702 (MYB) transcript levels, at 6 h to sequence miRNAs and mRNAs, and at 36 h to analyze cold tolerance.

The maximum photochemical efficiency of PSII (F_v/F_m) was measured with Portable Chlorophyll Fluorometer (PAM2500; Heinz Walz, Effeltrich, Germany), after the whole plants were dark adapted for 30 min. Minimal fluorescence (F_o) was measured during the weak measuring pulses and maximal fluorescence (F_m) was measured by a 0.8-s pulse light at 4,000 μmol m^-2 s^-1. F_v/F_m was determined using the third leaf of watermelon plants counting from the bottom to up. The calculation of F_v/F_m was done by the following formula as described by van Kooten and Snel (1990). F_v/F_m = (F_m - F_o)/F_m.

REL in the leaves was determined as previously described elsewhere (Zhou and Leul, 1998).

Eight independent small RNA libraries from plant leaves for four treatments (Control, CK; Melatonin, MT; Cold; Melatonin+Cold, MT-C), with two biological replicates for each treatment, were sequenced. The sequences have been deposited into the NCBI Sequence Read Archive database (SRP078211, SRA438995). Library construction and sequencing of small RNA were performed using a service provider Gene Denovo Co. (Guangzhou, China). Total RNA was extracted using miRNeasy Mini Kit (Cat#217004, QIAGEN GmBH, Germany) according to the manufacturer’s protocol. Total RNA integrity was measured on an Agilent 2100 Bioanalyzer system (Agilent) for quality control. A total of 16–35 nt RNA fragments were excised, purified from a PAGE gel, and ligated with 5′ and 3′ adaptors using T4 RNA ligase. Reverse transcription followed by PCR was used to create cDNA constructs based on the small RNA ligated with 3′ and 5′ adapters. Subsequently, the amplified cDNA constructs were purified from agarose gel, in preparation for sequencing analysis using the Illumina Genome Analyzer (Illumina, CA, USA) according to the manufacturer’s instructions.

The raw sequences were firstly processed by Illumina’s Genome Analyzer Pipeline software to filter out the adapter sequences, low quality as well as low-copy sequences. Then, the extracted small RNA sequences with 15–26 nt in length were subjected to cucurbit species mRNAs, Rfam¹ and Repbase² to discard mRNA, rRNA, tRNA, snRNA, snoRNA and repeat sequences. Finally, the remaining unique sequences were analyzed by BLAST against miRBase (Release 21)³. Solexa sequences with identical or related (one mismatch) sequences from mature miRNAs were identified as known miRNAs (Meyers et al., 2008).

To identify potential novel miRNAs in watermelon, rest of the unmapped small RNA sequences were searched by BLAST against watermelon genome downloaded from cucumber genome database⁴. The mappable sequences were then folded into a secondary structure using UNAfold software⁵. Only the non-coding sequences which could form a perfect stem-loop structure and meet the criteria for miRNAs prediction (Meyers et al., 2008) were then considered to be a potential novel miRNA candidate.

We chose the miRNAs increased or decreased by more than twofold and P-value <0.01 in two treatments as the criterion for a melatonin or cold response. Then, the heat map of differentially expressed miRNAs expression profile was drawn with MultiExperiment Viewer version 4.0 and clustering analysis was performed using a hierarchical clustering method (Eisen et al., 1998).

Eight independent mRNA libraries from plant leaves for four treatments (Control, CK; Melatonin, MT; Cold; Melatonin+Cold, MT-C), with two biological replicates for each treatment, were sequenced. The sequences have been deposited into the NCBI Sequence Read Archive database (SRP078211, SRA438977). Library construction and sequencing of mRNA were performed using a service provider Gene Denovo Co. (Guangzhou, China). After removal of low quality sequences, the clean reads were mapped to the watermelon reference genome, allowing up to one mismatch. The differentially expressed genes were identified using the R package edgeR (Robinson et al., 2010). The expression level of each unigene was calculated and normalized to generate FPKM (fragments per kilobase of exon per million mapped fragments). The false discovery rate (FDR) was used to determine the threshold of the P-value in multiple tests. In this study, the FDR < 0.05 and fold change ≥1.5 were used as significance cut-offs of the gene expression differences.

A plant miRNA target prediction server⁶, was used to predict putative miRNA target genes with default settings based on the library of Watermelon genome database, version 1. It reports all potential sequences complementary to an inquiring miRNA sequence with mismatches no more than a specified value for each mismatch type. The minimal weighed score ≤3.0 was applied in the prediction according to scoring schema of miRU by Zhang (2005). The other default settings as follows: length for complementarity scoring (hspsize), 20 bp; target accessibility-allowed maximum energy to unpair the target site (UPE), 25.0; flanking length around target site for target accessibility analysis, 17 bp in upstream/13 bp in downstream; range of central mismatch leading to translational inhibition, 9–11 nt (Brodersen et al., 2008).

Two cold-sensitive genes Cla020078 (CBF1) and Cla020702 (MYB) were chose to perform quantitative real-time PCR (qRT-PCR). Total RNA was extracted using a RNA extraction kit (Tiangen, Beijing, China) according to the supplier’s instructions. DNA contamination was removed using a purifying column. One microgram of total RNA was reverse-transcribed using the ReverTra Ace qPCR RT Kit (Toyobo, Osaka, Japan) following the supplier’s instructions. The gene-specific primers for qRT-PCR were designed based on their cDNA sequences, as follows: Cla020078 (F, AGCAGAGCCCTAACACAGGT; R, AATGGTCTTGAGTTGGG), Cla020702 (F, GATCCATTGACGGCACTAAC; R, TCGCTACAACGTCCTTCATC), and watermelon β-actin gene (F, CCATGTATGTTGCCATCCAG; R, GGATAGCATGGGGTAGAGCA) was used as an internal control (Kong et al., 2014). The qRT-PCR assays were performed using an iCycler Iq Multicolor PCR Detection System (Bio-Rad, Hercules, CA, USA). PCRs were performed using the SYBR Premix ExTaq II (2×) Kit (Takara, Tokyo, Japan). The PCR conditions consisted of denaturation at 95°C for 3 min, followed by 40 cycles of denaturation at 95°C for 30 s, annealing at 58°C for 30 s, and extension at 72°C for 30 s. The quantification of mRNA levels was based on the method of Livak and Schmittgen (2001).

The experiment was a completely randomized design with three replicates. Each replicate contained at least 10 plants. Analysis of variance was used to test for significance, and significant differences (P < 0.05) between treatments were determined using Tukey’s test.

Melatonin plays important regulatory roles in plant defense against various biotic and abiotic stresses (Zhang et al., 2014). In this study, we first analyzed the effects of exogenous melatonin on watermelon tolerance to cold stress. As shown in Figure 1, melatonin had a minimal effect on the leaf phenotypes, the maximum quantum yield of PSII (F_v/F_m), and REL at optimal growth temperatures. After 4°C treatment for 36 h, the blade edges of watermelons were wilted, and the F_v/F_m and REL were significantly reduced and increased, respectively. However, pretreatment with melatonin obviously alleviated cold-induced wilting of blade edges and the reduction and increase in leaf F_v/F_m and REL, respectively. For example, the F_v/F_m and REL decreased and increased by 15.45 and 19.82%, respectively in melatonin-pretreated plants after cold treatment, far less than the changes of 31.06 and 50.51% in control plants (Figure 1). Taken together, exogenous melatonin significantly enhanced watermelon tolerance to cold stress.

To analyze the effects of melatonin on the expression of cold-responsive genes and to choose the appropriate time point for miRNA expression profiling, we examined the dynamic responses of Cla020078 (CBF1) and Cla020702 (MYB) to melatonin or/and cold treatment. Transcript levels of CBF1 and MYB were slightly induced by melatonin at optimal growth temperatures (Figure 2). After cold treatment (4°C), they were rapidly upregulated, reaching peak levels at 6 h, and subsequently declining to original levels at 36 h. Furthermore, pretreatment with melatonin improved the induction of CBF1 and MYB expression by cold stress. For instance, at 6 h after cold treatment, the expression levels of CBF1 and MYB in melatonin-pretreated plants were upregulated by 5.99- and 85.69-fold, respectively, far more than the increases of 3.18- and 43.38-fold observed in control plants.

To examine whether miRNAs are involved in melatonin-mediated cold response in watermelon, we performed high-throughput sequencing analysis of miRNA in watermelon leaves treated with and without the melatonin and cold stress. A total of 9,704,684/10,831,539 (Control, CK-1/2), 10,124,858/12,383,237 (Melatonin, MT-1/2), 12,812,180/9,718,613 (Cold-1/2), and 9,806,063/11,489,173 (Melatonin+Cold, MT-C-1/2) raw reads were obtained (Supplementary Table S1). After the removal of rRNAs, tRNAs, snRNAs, and snoRNAs, a total of 5,327,079/5,914,071, 5,635,109/6,581,846, 6,219,825/5,085,935, and 5,409,138/6,272,887 sRNA sequences for CK-1/2, MT-1/2, Cold-1/2, and MT-C-1/2 were obtained, respectively. The majority of redundant reads were in the range of 19–26 nt, among which the most abundant sequences were 24 nt long in all libraries (Figure 3A). To identify miRNAs in watermelon, all sRNA sequences were compared to known plant miRNAs in miRBase Release 21. A total of 440 known unique miRNAs with high sequence similarity to known plant miRNAs were identified (Supplementary Table S2). By mapping all unique sRNA sequences to the watermelon genome and predicting the hairpin structures for their flanking sequences, 106 novel miRNAs candidates were identified. The expression level of each miRNA was normalized to generate TPM (tags per million) and the correlation of detected miRNA expression between two biological repeats for each treatment was analyzed. As shown in Figure 3B, there was a strong correlation (R = 0.99, P < 0.0001) between the two biological replicates for all four treatments, suggesting that the sequencing results are highly reliable.

A total of 50 miRNAs were significantly up- or downregulated by melatonin under unstressed or cold conditions (Figure 4). Among them, four miRNAs such as miR170-5p, novel-m0074-3p, novel-m0044-3p, and novel-m0040-3p were upregulated, while one miRNA (miR399-5p) were up- and downregulated by melatonin pretreatment, under both unstressed and cold conditions. Ten and nine miRNAs were specifically up- and downregulated by melatonin at optimal growth temperatures, respectively. Twelve and 20 miRNAs were specifically up- and downregulated by melatonin after cold treatment for 6 h, respectively. A total of 54 miRNAs were significantly up- or downregulated by cold stress in control or melatonin-pretreated plants, respectively. Among them, four miRNAs (miR399-5p, miR170-3p, miR482-5p, and novel-m0031-5p) and three miRNAs (miR159-5p, novel-m0033-5p, and novel-m0038-5p) were up- and downregulated by cold stress, respectively, in both control and melatonin-pretreated plants. In total, 18 and 12 miRNAs were specifically up- and downregulated by cold in control plants, respectively. While 8 and 12 miRNAs were specifically up- and downregulated by cold in melatonin-pretreated plants, respectively.

To examine the biological functions of differentially expressed miRNAs, the miRNA sequences were searched against watermelon genomic sequences using the plant miRNA potential target finder⁷ to predict target mRNAs. miRNAs negatively regulate mRNAs by inducing their cleavage or repressing their translation (Bartel, 2004). In this study, a total of 505 mRNAs were predicted to be cleaved by miRNAs. Then, we performed RNA-seq analysis in watermelon leaves treated with and without the melatonin and cold stress using the high-throughput Illumina Solexa system. Reads from each sample were aligned to the C. lanatus reference genome. The expressions of predicted mRNAs in comparison of MT/CK, MT-C/Cold, Cold/CK, and MT-C/MT were analyzed. The expression patterns of only 49 detected mRNAs were opposite with respect to those of 21 corresponding miRNAs (Table 1). The negative correlation coefficients between miRNAs and their target mRNAs regulated by melatonin and cold were -0.74 (P < 0.001) and -0.48 (P < 0.001), respectively, suggesting that these 49 target mRNAs were likely cleaved by the corresponding miRNAs (Figure 5). Additionally, 34 miRNAs were searched against a total of 87 potential target mRNAs via translational inhibition (Supplementary Table S3). The functions of these mRNAs were annotated in the Cucurbit Genomics Database⁸.

Some of the 27 potential target mRNAs cleaved by melatonin-responsive miRNAs were induced or repressed by melatonin at optimal growth temperatures, but not under cold stress (Figure 6A; Supplementary Table S4). Interestingly, most miRNAs including some cold-responsive genes, such as Cla011198 (BHLH), Cla017355 (WRKY), and Cla022458 (DREB, dehydration-responsive element-binding gene), were upregulated by melatonin during cold stress, but were not affected by melatonin under unstressed conditions. Most of 25 potential target mRNAs cleaved by cold-responsive miRNAs were repressed or unaffected by cold stress, except for Cla017745 (LEA) and Cla017373 (R2R3-MYB), those were upregulated (Figure 6B). However, melatonin pretreatment alleviated cold-induced repression of target mRNAs or improved cold-induced increases in target mRNAs. For instance, LEA and R2R3-MYB in melatonin-pretreated plants were upregulated by 5.53- and 14.91-fold under cold stress, far more than the changes of 2.82- and 2.44-fold in control plants, respectively. MYB, MDAR (monodehydroascorbate reductase), and PPRP (pentatricopeptide repeat protein) were significantly repressed by cold stress in control plants, but were minimally affected in melatonin-pretreated plants. Accordingly, a set of potential target mRNAs cleaved by corresponding miRNAs might be involved in the response of melatonin-pretreated plants to cold stress.

We also subjected the predicted target mRNAs to Gene Ontology classification based on their involvement in process in Cucurbit Genomics Database⁹ with watermelon 97103 v1 (Figure 7; Supplementary Figure S1). The target genes inhibited by miRNAs via cleavage and translational repression were involved in various biological processes, with 19 and 38 mRNAs involved in cellular process regulation, 18 and 31 mRNAs involved in responses to biotic or abiotic stress, and 5 and 10 mRNAs involved in signal transduction, respectively. Thus, a number of potential target mRNAs might be involved in watermelon responses to cold stress by direct or indirect ways.

In recent years, melatonin has emerged as a focus of research in plant science. It functions in many aspects of plant growth and development, and regulates plant tolerance to abiotic stresses, such as salinity, drought, radiation, excess copper, and chemical stresses (Zhang et al., 2014). Recently, several studies have shown that exogenous melatonin with optimum concentrations enhances cold tolerance in Arabidopsis and Triticum aestivum L. (Bajwa et al., 2014; Shi and Chan, 2014; Turk et al., 2014), since the initial observation that exogenous melatonin attenuates cold-induced apoptosis in carrot suspension cells (Lei et al., 2004). In the current study, cold stress caused wilted blade edges, decreased maximum photochemical efficiency of PSII (F_v/F_m), and increased REL, which indicates membrane damage (Figure 1). Notably, in agreement with previous studies, the adverse effects of cold were significantly alleviated by the application of melatonin. Moreover, melatonin treatment up-regulated two cold-responsive genes Cla020078 (CBF1) and Cla020702 (MYB) in watermelon, which is in line with previous studies in Arabidopsis (Bajwa et al., 2014; Shi and Chan, 2014), suggesting that melatonin activates signaling pathways during cold stress via regulation of cold-responsive genes (Figure 2).

Transcriptional control of the expression of cold-responsive genes is well known (Chinnusamy et al., 2010), but miRNAs, which could regulate stress responses by controlling the expression of cognate target genes, have recently been added to the suite of cold-responsive gene regulatory networks. Comparative profiles of miRNA expression among various plant species (Arabidopsis, Brachypodium, and Populus) during cold stress revealed similarities and differences in miRNA regulation (Liu et al., 2008; Lu et al., 2008; Zhang et al., 2009). In this study, 440 known miRNAs and 106 novel miRNAs were identified using Solexa high-throughput sequencing (Supplementary Table S2). Based on an analysis of differentially expressed miRNAs, we identified a set of miRNAs including 16 known miRNAs and 21 novel miRNAs that were significantly up- or downregulated by cold stress, which exhibited similarities and differences with other species (Figure 4; Khraiwesh et al., 2012; Sunkar et al., 2012). The discrepancies in cold-responsive miRNA expression might be attributed to the differential responses of unrelated plant species and variation in experimental conditions and methods. Importantly, pretreatment with melatonin affected the expression of a set of miRNAs under unstressed or cold conditions, and changed the expression of cold-responsive miRNAs. These results suggest that a number of miRNAs are involved in melatonin-mediated enhancement in cold tolerance in watermelon. To the best of our knowledge, this is the first report on the interaction between melatonin and miRNAs in regulation of plant tolerance to stress, although a few studies have reported such interactions in animal cell lines (Lee et al., 2011; Sohn et al., 2015).

To examine the functions of melatonin- or cold-responsive miRNAs, we predicted their target mRNAs. We detected 49 mRNAs with opposite expression patterns with respect to their corresponding miRNAs, and inferred that these mRNAs are likely cleaved by the corresponding miRNAs (Table 1). Among the 49 potential target mRNAs, several mRNAs were closely associated with plant defense against cold stress by activating signaling transduction or encoding functional proteins. These mRNAs include Cla001826, Cla011198, Cla022458, Cla004188/Cla017373, Cla017355, Cla017745, and Cla018379 which encode CDPK, BHLH transcription factor, DREB protein, MYB transcription factor, WRKY transcription factor, LEA protein, and MDAR, respectively. As Ca²⁺ sensors, CDPKs those can directly bind Ca²⁺, are major players in coupling cold stress signals to specific protein phosphorylation cascades, leading to the activation of downstream cold-responsive genes, such as ICE1, a MYC-type BHLH transcription factor (Janská et al., 2010). ICE1 is activated by low temperatures and triggers CBF/DREB expression, which promotes the accumulation of cold-responsive gene products such as LEA, antioxidant enzymes, and molecular chaperones (Janská et al., 2010). MYB and WRKY transcription factor gene families have been suggested to play important roles in the regulation of transcriptional reprogramming associated with plant responses to cold stress (Shinozaki et al., 2003; Zhu et al., 2005; Chen et al., 2012). Overexpression of some MYB or WRKY genes enhances plant tolerance to cold stress and activates the expression of a set of cold-responsive genes such as LEA, DREB2A, and genes encoding antioxidant enzymes (Dai et al., 2007; Niu et al., 2012; Yang et al., 2012). LEA is important for membrane stabilization and the protection of proteins from denaturation when the cytoplasm becomes dehydrated (Janská et al., 2010). MDAR is an important antioxidant enzyme that is critical to alleviate cold-induced oxidative damage (Chen and Paul, 2002). Importantly, these defense-related genes were significantly induced by melatonin, exhibiting contrasting expression patterns with respect to the corresponding miRNAs, including miR159-5p, miR858, miR8029-3p, and novel-m0048-3p (Figures 5 and 6). Accordingly, it could be inferred that some cold resistance pathways in melatonin-pretreated plants were promoted by a number of downregulated miRNAs under cold stress.

Some other potential target genes might also be involved in watermelon responses to cold stress by direct or indirect ways. For example, low temperatures can inhibit the uptake and transport of nitrate (Warren, 2009; Leffler et al., 2013; Dechatiwongse et al., 2014). The significantly upregulated expression of Cla015456 (NRT) by melatonin during cold stress suggests that the active transport of nitrate could be helpful in the acclimation of plants under adverse environmental conditions (Zhang et al., 2005). Additionally, the differentially expressed miRNAs were predicted to target a set of mRNAs via translational repression (Supplementary Table S3). Some of these mRNAs, such as MYB, MYC, and CYP450, are important in transcriptional regulation and protection in cold responses, respectively. However, whether these predicted target mRNAs are regulated by corresponding miRNAs via translational inhibition requires further studies.