Comparative transcriptomics reveals new insights into melatonin-enhanced drought tolerance in naked oat seedlings

Plant materials and treatments

The experiment was conducted in a smart greenhouse at Zhangjiakou Academy of Agricultural Sciences, Zhangjiakou, Hebei Province, China. This experiment utilized two naked oat varieties: ‘Huazao No.2’ (H2) and ‘Jizhangyou No.15’ (J15). The Zhangjiakou Academy of Agricultural Sciences developed both H2 and J15. Naked oat seeds were germinated in an incubator at 25 °C in darkness for 24 h. The germinated seeds were sowed in pots (30 cm long, 18 cm wide, 12 cm high) filled with a 3:1:1 (v:v:v) mixture of topsoil (sampled from the 0–20 cm upper soil layer; organic matter content 10.4 g kg⁻¹, total N 1.31 g⁻¹, alkali-hydrolyzable N 41.6 mg kg⁻¹, available phosphorus 11.5 mg kg⁻¹, available potassium 121 mg kg⁻¹), nutrient soil (Pindstrup Plus, Ryomgård, Denmark; pH 6.0, screened to 0–6 mm), and sand in an environmentally controlled greenhouse (day/night air temperature, 25/22 °C; photoperiod, 14/12 h; 600 µmol photons m⁻² s⁻¹ light during the day; relative humidity, 60–70%).

Plants of each variety were divided into two groups based on melatonin (MT) application and soil relative water content (SRWC) as follows: well-watered (CK, set as the control group), 75 ± 5% SRWC; drought stress (DS), 45 ± 5% SRWC Xiao et al. (2020); drought stress with melatonin (DS+MT), 45 ± 5% SRWC and sprayed with 100 µM MT (Sigma-Aldrich, MO, USA).

Ten days after emergence, plants were subjected to drought treatments and MT applications. The method of spraying was as follows: all leaves were wet but not dripping; all other treatments received sprayed distilled water; plants were sprayed five times, with each spray given every other day. The SRWC was controlled by the weighing method, and plants were weighed every other day and watered to maintain the appropriate SRWC range throughout the experimentation. Each experiment had nine replicates.

Measurement of plant growth and soil and plant analyzer development values

The phenotypes of naked oat seedlings were determined starting 5 days after drought treatment (DAD). The aboveground phenotype of naked oats was determined every 10 days, three times. Plant height was measured using a ruler. Leaf area was calculated by the length-width coefficient method. The chlorophyll concentrations of the leaves were determined with a soil and plant analyzer development (SPAD) meter (SPAD-502, Konica-Minolta, Tokyo, Japan).

Measurement of leaf relative water content

Leaf relative water content (LRWC) in functional leaves was measured at 25 DAD, according to the method described by Xiao et al. (2020). Briefly, the fresh weight (FW) measures were immediately made after sampling. Then, samples were immersed in distilled water for 4 h at room temperature (25 °C) before measuring their turgid weight (TW). The leaf samples were then blotted dry and weighed after being oven-dried at 85 °C for 48 h, after which their dry weight (DW) was measured. The LRWC was calculated based on the following formula: LRWC (%) = [(FW − DW)/(TW − DW)] × 100%.

Gas exchange and instantaneous water use efficiency

Three representative oat seedlings were selected from each plot for measuring net photosynthetic rate (Pn) and transpiration rates (Tr) using an LI-6400 portable photosynthesis system (LI-COR, NE, USA) under ambient conditions (irradiation = 600 µmol m⁻² s⁻¹, leaf temperature = 25 ± 1 °C) at 25 DAD. Instantaneous water use efficiency (WUE_inst) was calculated from the above indicators using the formula: WUE_inst = Pn/Tr.

RNA extraction and cDNA library synthesis

Leaf samples from CK, DS, and DS +MT, were collected at 25 DAD and immediately stored at −80 °C. Each treatment had three biological replicates and no technical replicates. The labeled frozen tissues were ground to a fine powder in liquid N₂, and the total RNA was isolated using the TRIzol reagent (Invitrogen, Waltham, MA, USA). The RNA quality was evaluated on a 1% (w/v) agarose gel electrophoresis and the RNA Nano 6000 Assay Kit of the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The mRNA was extracted from total RNA using poly-T oligo-attached magnetic beads. The mRNA was fragmented using divalent cations under elevated temperature in the First-Strand Synthesis Reaction Buffer (5X). The library fragments were purified with the AMPure XP system (Beckman Coulter, Brea, CA, USA), preferentially selecting 370∼420 bp fragments. The selected fragments were amplified by PCR and purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA) to obtain the final library. The cDNA library was quantified using a Qubit2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and diluted to 1.5 ng/ul. The library insert sizes were determined using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). After extracting the inserts with preferred sizes (370∼420 bp), a qRT-PCR analysis accurately quantified the library concentration (>2 nM) to ensure quality sequencing of the libraries. Then, high-quality cDNA libraries were pair-end sequenced on the Illumina HiSeq 6000 platform (Illumina, San Diego, CA, USA). The transcriptome datasets are available at the NCBI Sequence Read Archive (SRA), accession numbers SRX12257066–SRX12257083 (Bioproject number PRJNA764537).

Sequencing read processing and mapping

The high-throughput image data were automatically converted into sequence data (reads) using the CASAVA base recognition software (Illumina, San Diego, CA, USA). Raw reads were cleaned using fastp (version = 0.20.1) with default parameters. First, adapter sequences were trimmed off from all raw reads. Next, undetermined bases (N) were removed, and low-quality reads (Phred score ≤ 20 covering > 50% of the whole read length) were dropped. The Q20, Q30, and GC contents of clean data were calculated by (version = 2.2.1), and clean reads were retained for downstream analysis. After that, the clean reads were de novo assembled using the Trinity software (v2.6.6) (Garber et al., 2011). The BUSCO software evaluated the splicing quality, accuracy, and completeness of the assemblies from Trinity.fasta, unigene.fa, and cluster.fasta according to the assembly proportion and completeness. Mapping was performed using the FeatureCounts v1.5.0-p3. Filtered Gene Sets were identified, and the resultant read counts were normalized by FPKM (Fragments Per Kilobase of transcript per Million mapped reads) to measure the transcript expression. Gene function was annotated based on the following databases: NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences(Nt), and UniProt (Swiss-Prot) with a 10⁻⁵e-value threshold.

Differential expression analysis

Differentially expressed genes (DEGs) between CK, DS, and DS+MT were identified using the DESeq2 R package (1.20.0) via the Novomagic v3.0 platform (https://magic.novogene.com), with a false discovery rate (p_adj) < 0.05 and log₂(ratio) ≥ 1 (Li et al., 2021) thresholds for significant differential expression. Enrichment of Gene Ontology (GO) terms was determined using agriGO (https://bio.tools/agrigo), and enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was determined using KOBAS 3.0 (http://kobas.cbi.pku.edu.cn/kobas3/?t=1). The significantly enriched GO items and KEGG pathways with FDR-corrected p-value ≤ 0.05 were analyzed and visualized using OmicShare tools, a free online platform for data analysis (https://www.omicshare.com/tools).

Quantitative real-time PCR (qRT-PCR) validation

Nine DEGs obtained from the Illumina RNA-seq were randomly selected for qRT-PCR validation with Actin as the internal reference gene. Primer3Plus (https://www.primer3plus.com/) was used to create gene-specific primers (Table S1). Each sample had three biological and technical replicates. The quantitative real-time PCR (qRT-PCR) was performed on an ABI 7500 Real-Time PCR System (Applied Biosystems, CA, USA) following the manufacturer’s instructions. The amplification program involved: 95 °C for 2 min followed by 45 cycles of 95 °C for 5 s and 60 °C for 1 min. The relative expression level of each unigene was calculated by the 2^−ΔΔCT method (Livak & Schmittgen, 2001). The 7500-system SDS Dissociation Curve Analysis Software examined the amplification specificity of each run.

Statistical analyses

Statistical analyses for transcriptomic data were performed using three biological replicates and five biological replicates for the other analyses. Analysis of variance (ANOVA) was performed using the SPSS 26.0 software (IBM Corp., Armonk, NY, USA) following Duncan’s multiple range tests at the 5% probability level. The results were obtained using SPSS and GraphPad software and presented as mean ± standard deviation (SD) of three independent biological experiments. The RNASeqPower package in R generated >0.91 size estimates for a 2.0 fold change, verifying the high-quality data sets used for differential gene expression analysis.

Shoot characteristics under exogenous melatonin

Drought stress significantly affected the growth of the two varieties of naked oat seedlings (Figs. 1A and 1B). Spray application of MT alleviated the effects of drought stress compared with the DS treatment at 25 DAD, especially for H2, which had more upright leaves. In H2, plant height (PH), leaf area (LA), dry weight (DW), and fresh weight (FW) of plants under drought stress decreased by 18.85, 52.35, 43.72, and 57.17%, respectively (P < 0.05). However, J15 registered 13.35, 27.03, 32.59, and 45.19% decrease (P < 0.05), respectively (Figs. 1C, 1D, 1E, 1E and 1F). The shoot characteristic data show that MT spraying significantly alleviated drought stress in H2, as PH, LA, DW, and FW increased by 9.48, 66.15, 53.42, and 47.6%, respectively, than DS treatments (P < 0.05). Similarly, spraying MT improved the shoot characteristics of J15 plants.

Exogenous melatonin enhances photosynthetic indexes and water use efficiency

Drought stress significantly decreased the Pn and SPAD value in DS-treated leaves (P < 0.05) (Figs. 2A and 2B). In H2 plants, the SPAD value and Pn were significantly greater in melatonin-treated than DS treated leaves (P < 0.05). However, the MT mitigation effect on J15 was nonsignificant. The LRWC reflection of each treatment was consistent with the SPAD and Pn results (Fig. 2C). The LRWC reflection of each treatment was consistent with the SPAD and Pn data, contrary to WUE_inst. Drought stress increased the WUE_inst of H2 and J15 leaves by 77.03 and 28.16% than CK treatment (P < 0.05) (Fig. 2D). MT treatment further increased WUE_inst, but not significantly.

Global description of naked oat transcriptome

Eighteen cDNA libraries were prepared from seedlings under CK, DS, and DS + MT treatments to unveil the potential regulatory mechanisms underlying melatonin-mediated drought responses in naked oats. A total of 22.40–27.21 million raw reads were obtained from each cDNA library. After filtering the raw data and checking both the sequencing error rate and GC content distribution, 274.68 Gb of clean reads were obtained, with an average of 15.26 Gb of reads per sample. Per sample, Q20 and Q30 were 97.85 and 93.38%, respectively (Table S2).

Cleaned reads were de novo assembled with Trinity, yielding 302,169 transcripts and 95,202 unigenes with 2,110 and 1,239 bp N50 lengths, 726 and 510 bp N90 lengths, respectively (Table S3).

The length distribution of assembled oat transcripts and unigenes revealed that 25.30 and 18.68% of the total transcripts and unigenes were greater than 2,000 bp, respectively (Fig. S1).

All assembled oat unigenes were validated and annotated by implementing similarity searches against several public databases (Fig. S2). Of the 95,202 unigenes, 46.15 and 43.55% aligned to the NR and NT databases. Moreover, 14.65, 35.27, and 7.59% of unigenes exhibited significant matches with the KO, GO, and KOG databases. Ultimately, 57,354 unigenes (60.24%) were annotated in at least one database.

We analyzed the species distribution of the All-Unigene datasets by aligning the sequences with the NR database. A species distribution map was drawn based on the results to characterize the sequence similarity between naked oat and other species (Fig. S3). Naked oat genes were most similar to Aegilops tauschii (32.9%). Over 79.3% of the distinct All-Unigenes sequences had top matches with sequences from other Gramineae plant species. Pearson correlation analysis indicated that the three biological replicates had highly consistent transcriptome profiles across all tissue types (Fig. S4).

Differential gene expression profiling

We compared the expression profiles of both varieties under CK, DS, and DS+MT treatments at 20 DAD to obtain detailed information about the expression profiles of genes. Uniquely aligned reads were used to estimate gene expression levels as FPKM. As illustrated by the Venn diagram, only significant DEGs were considered for analyzing expression patterns (Fig. 3).

The total DEGs between DS and CK treatments was 4,672 for H2, with 3,586 up-regulated and 1,086 down-regulated (Fig. 3A). For H2 leaves, 4,672 genes were differentially expressed between DS and CK treatments. There were 1,984 DEGs between DS+MT and CK treatments, half the number of DEGs between DS and CK treatments. Between the DS+MT and DS treatments, there were 3,418 DEGs, of which 1,635 were up-regulated and 1,783 were down-regulated. In J15 leaves, 1,160 genes were differentially expressed between DS and CK treatment, 2,661 between DS+MT and CK treatments, and 4,654 between DS+MT and DS, of which 2,073 were up-regulated 2,581 were down-regulated (Fig. 3B). Thus, different varieties responded differently to drought and melatonin treatments. As shown in the Venn diagram (Figs. 3C and 3D), 342 up-regulated and 123 down-regulated genes were observed in H2 and J15 varieties between the DS +MT and CK treatments.

The DEGs between DS+MT and DS are presented as the number of co-up-regulated (39) or co-down-regulated (89) genes in both varieties (H2 and J15), respectively (Fig. 3E). Interestingly, there are 172 genes with opposite trends in both varieties, where 41 genes were up-regulated in H2 and down-regulated in J15. Besides, 131 genes were down-regulated in H2 and up-regulated in J15. GO and KEGG analyses were performed for the DEGs between DS+MT and DS.

Screening and GO classification of melatonin-induced DEGs

GO enrichment analysis was performed on the transcriptome data to understand the role of melatonin supplementation in naked oat seedlings under drought stress. GO terms classified the enrichment of genes into biological processes, molecular function, and cellular components (Fig. 4).

In this research, we focused on the results of two naked oat varieties under drought stress with supplemental melatonin. The results obtained from the GO enrichment analysis with the genes included in each category are summarized in Table S4. Under biological process, the enriched terms included “response to water”, “response to acid chemical”, “response to oxygen-containing compound”, and “response to the abiotic stimulus”. For molecular function, the enriched terms included “glucosidase activity” and “oxidoreductase activity”. The cellular component included “integral component of Golgi membrane”, “intrinsic component of Golgi membrane”, “glucosidase II complex”, and “glycine reductase complex”. These results suggest that leaves of drought-stressed naked oat seedlings underwent global transcriptional reprogramming under supplementation with melatonin.

KEGG pathway enrichment analysis

DEGs were subjected to KEGG enrichment analysis to investigate which pathways were involved in melatonin supplementation under drought stress in naked oats. We identified five statistically significant enriched pathways, including “Plant hormone signal transduction”, “beta-Alanine metabolism”, “Nitrogen metabolism”, “alpha-Linolenic acid metabolism”, and “Linoleic acid metabolism” (q-value < 0.05) (Fig. 5, Table S5).

The “plant hormone signal transduction pathway” (PHSTP) contained the highest number of DEGs, including 16 up-regulated and two down-regulated genes. Among the 16 up-regulated genes, 15 were concentrated in region III of the Venn diagram, and the two down-regulated genes belonged to regions II and III, respectively (Fig. 3). Seventeen DEGs enriched in this pathway were related to abscisic acid (ABA), and one was related to indole acetic acid (IAA) (Table S6). These DEGs observed in the DS + MT treatment were generally expressed at lower levels than under the DS treatment (Table 1). Four valuable down-regulated genes were concentrated in the nitrogen metabolism pathway, including the nitrate reductase, nitrate transporter, and Alpha carbonic anhydrase genes (Table S7). These DEGs observed in the DS + MT treatment were generally expressed at lower levels than under DS treatment.

Expression verification of naked oat seedling DEGs

The reliability of RNA-Seq gene expression was confirmed by qRT-PCR. Nine DEGs consisting of six up-regulated and three down-regulated genes were selected for qRT-PCR analysis (Fig. 6). The fold changes obtained from RNA-seq and qRT-PCR results were highly correlated (r = 0.998, p < 0.01) (Fig. S5). Hence, the present study’s naked oat seedling transcriptome analysis was reliable.