Transcriptomic and morphophysiological responses of Selenicereus undatus to single and combined drought and nickel stress: melatonin-mediated mitigation

Selenicereus undatus (2n = 2x = 22), commonly known as pitaya or dragon fruit, is a tropical epiphytic cactus native to Mexico that has become an important specialty crop in Hainan and is extensively cultivated across South China, Southeast Asia, and Latin America^[1]. Valued for its nutritional profile and economic significance, pitaya is rich in antioxidants and essential minerals while being low in calories, fat, and protein^[2]. However, abiotic stresses threaten pitaya yields by impairing growth and development at both the seedling and maturity stages. Although the crop activates various defensive strategies to mitigate stress-induced injury^[3], severe conditions can still cause substantial yield losses^[2].

Heavy metal contamination in agricultural soils has increased significantly over the past two decades, with marked exceedance rates observed since 2017 due to rapid urbanization and industrial runoff^[4]. Nickel (Ni) has become an escalating concern in agricultural systems worldwide, particularly in tropical and subtropical regions. The primary sources of Ni contamination are anthropogenic, including industrial discharges, mining activities, and intensive agricultural practices, with natural sources contributing minimally^[5]. Plants primarily absorb Ni through roots and, to a lesser extent, leaves. Once absorbed, Ni can enter the food chain, creating serious health risks, especially to vulnerable populations^[6]. In Sanya soil, Ni is 6.56-fold higher than the background value (7.24−27.69 mg·kg^–1)^[7]. Trace Ni (< 500 µg·kg^–1) supports plant growth but is toxic at higher levels^[8]. Although plants sequester heavy metals in vacuoles and cell walls, excess Ni continues to impair photosynthesis and reduce biomass, resulting in significant yield losses^[9,10]. As such, enhancing the tolerance of economically important fruit crops, such as pitaya, to Ni stress has emerged as a critical research priority.

Similarly, changing climatic conditions have led to widespread environmental issues, with rising temperatures intensifying droughts worldwide^[11]. In parallel, warming-driven drought intensifies reactive oxygen species (ROS) formation, compromising photosynthesis and membrane integrity; consequently, antioxidant systems become pivotal. In response to drought, pitaya accumulates ROS, which damage cellular components and metabolic pathways, thereby impairing growth and development^[12]. The ROS scavenging mechanisms associated with antioxidants, such as superoxide dismutase (SOD) and catalase (CAT), have been shown to protect cell membranes from oxidative damage caused by osmotic stress. Stomatal regulation under water deficit can alter transpiration-driven fluxes and thereby influence ion transport and metal uptake. In pitaya, which exhibits crassulacean acid metabolism (CAM)-like physiology, stomatal behavior and carbon assimilation under water deficit may differ from those of typical C3 crops; therefore, drought responses should be interpreted in the context of this water-conserving photosynthetic habit^[10]. The plant hormone abscisic acid (ABA) plays a crucial role in enhancing drought tolerance, primarily by regulating guard cells to keep stomata closed and maintain water retention^[13]. Melatonin has recently gained attention as a promising technological intervention for mitigating heavy metal and drought stress in plants^[12,14]. Melatonin usually limits ROS accumulation, regulates ion homeostasis, delays senescence, and promotes root and bud growth^[15,16].

However, to the best of our knowledge, research on melatonin-mediated mitigation of nickel and drought stress in pitaya remains limited, and the underlying responses are still unclear. To address this gap, we evaluated the morphophysiological and transcriptomic responses of pitaya under single and combined drought and nickel stresses. We also examined the role of exogenous melatonin in alleviating these imposed stress conditions. In addition, we analyzed transcriptome profiles and co-expression gene networks using RNA-seq and weighted gene co-expression network analysis (WGCNA) to identify candidate genes and regulatory pathways associated with stress responses. The recent availability of chromosome-level genome assemblies of multiple pitaya species has also provided an improved genomic framework for transcriptome-based studies in the genus Selenicereus^[17−19]. This study provides new insights into melatonin-associated stress responses in pitaya and offers a practical approach for studying the stress adaptation of other Cactaceae plants.

The red-peel white-flesh pitaya (Selenicereus undatus) seeds were surface-sterilized and germinated in Petri dishes under controlled conditions at the Sanya Nanfan Research Institute, Yazhou Bay Seed Science & Technology City, Hainan University, China. Seedlings exhibiting uniform growth were then transplanted into half-strength Hoagland solutions^[20]. The entire cultivation process was conducted in a controlled growth chamber under a 16-h light/8-h dark period, a constant temperature of 24 °C, and 70% relative humidity. The treatments comprised eight groups: CK (control), M (melatonin, 100 µM), D (20% w/v PEG-6000), DM, Ni (NiCl₂, 125 µM), NiM, DNi, and DNiM. Melatonin (HPLC ≥ 98%, Cat. No. B21269, YuanYe Bio., Shanghai, China) was dissolved in ethanol to prepare a 1 mM stock solution, protected from light, and stored at 4 °C. A working solution was prepared by diluting the stock into the hydroponic nutrient solution to a final concentration of 100 µM melatonin. The selected concentrations were intended to impose clear but non-lethal stress while allowing evaluation of melatonin-associated mitigation. Specifically, 100 µM melatonin and 125 µM NiCl₂ were selected in line with previous pitaya stress studies^[2,21], whereas 20% PEG-6000 was used here as a reproducible osmotic treatment to impose drought-like stress under the hydroponic conditions of the present study. For morphological and biochemical assays, three biological replicates were used per treatment (24 plants across eight treatments). For RNA-seq, two biological replicates per treatment were collected (16 RNA-seq libraries) and then rapidly cryopreserved at –80 °C.

Phenotyping and Ni quantification

At the end of 60 d, three biological replicates per treatment were harvested for phenotypic analysis. Root length (RL), shoot length (SL), and cladode length (CL) were measured in cm using a ruler from the base to the tip of the respective organ. Cladode diameter (CD) was measured at the central region of the cladode using a vernier caliper in mm. Root weight (RW) and shoot weight (SW) were determined immediately after harvest using an analytical balance in grams (Supplementary Table S1). The Ni content in dried seedling tissue was determined by graphite furnace atomic absorption spectrophotometry (GFAAS; GTA-120) after nitric acid digestion^[22]. Results are expressed as µg Ni g^–1 dry weight.

Antioxidant enzyme assays and chlorophyll determination

Antioxidant enzyme activities and chlorophyll contents were determined using commercial assay kits from Nanjing Jianchen Bioengineering Institute (www.njjcbio.com/) following the manufacturer's protocols. Fresh cladode tissue was collected, immediately frozen in liquid nitrogen, and ground to a fine powder. For enzyme assays, 0.5 g of tissue was homogenized on ice in the extraction buffer specified by each kit, centrifuged to obtain the supernatant, and the supernatant was used for spectrophotometric measurements. Ascorbate peroxidase (APX; kit A123-1-1), catalase (CAT; Kit A007-1-1), and superoxide dismutase (SOD; Kit A001-4-1) were measured at 290, 405, and 550 nm, respectively, and reported on a fresh weight basis. Chlorophyll was determined using the chlorophyll assay kit (A147-1-1). Briefly, 0.1 g of fresh cladode tissue was extracted under light-protected conditions, adjusted to 5 mL with extraction reagent, incubated in the dark until the tissue residue became colorless (≥ 3 h), centrifuged, and the absorbance was measured at 645 and 663 nm to calculate chlorophyll a and chlorophyll b (mg·g^–1 FW).

RNA extraction and library preparation

Total RNA was extracted from pitaya cladode tissue using the RNAprep Pure Plant Kit (Tiangen Biotech, Beijing, China) following the manufacturer's protocol. RNA integrity and concentration were assessed before library construction, and 3 µg total RNA per sample was used as input for library construction. Poly(A)+ mRNA was enriched from total RNA using the MicroPoly(A) Purist Kit (Ambion, USA) according to the manufacturer's instructions. Strand-specific cDNA libraries were then constructed using the NEBNext Ultra RNA Library Prep Kit for Illumina, including RNA fragmentation, first- and second-strand cDNA synthesis, end repair and dA-tailing, adaptor ligation, size selection, and PCR enrichment. Library quality and concentration were evaluated before sequencing.

Sequencing, alignment, and quantification

Libraries were sequenced by Biomarker Technologies Co. Pvt. Ltd on the Illumina HiSeq™ 2500 platform to generate paired-end reads. Raw FASTQ files were processed using NGS QC Toolkits v2.3 to trim adapters, poly-N regions, and low-quality reads, yielding clean reads for downstream analyses. Clean reads were aligned to the reference genome^[17] with HISAT2 v2.0.5^[23]. Alignments were quantified at the gene level using StringTie to obtain the raw counts^[24]. DESeq2 v1.20.0 was used for differential expression under a negative binomial model. p-values were adjusted by Benjamini–Hochberg (BH) and reported as FDR. Differentially expressed genes (DEGs) were defined as FDR < 0.05 and |log₂FC| ≥ 3; the stringent fold-change cutoff was applied to prioritize robust, high-confidence expression changes. FPKM (fragments per kilobase per million) values were computed to visualize sample-to-sample Pearson correlation and to perform hierarchical clustering, and all DE inferences used raw counts. Functional annotation was performed with BLASTx against the Swiss-Prot protein database (EMBL-EBI, Cambridge, UK), the non-redundant NCBI nr protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa Laboratories, Kyoto, Japan), the Gene Ontology (GO) database, and the Cluster of Orthologous Genes (COG) database (NCBI, Bethesda, MD, USA). GO/KEGG enrichment was performed with ClusterProfiler v4.0 using the expressed gene set as background, with significance at FDR < 0.05.

Co-expression network analysis (WGCNA)

Co-expression network analysis was performed using the WGCNA v1.47 package in R^[25]. To reduce noise, genes with low expression (FPKM ≤ 2) were excluded. An unsigned, weighted network was constructed from FPKM expression matrices using a soft-thresholding power β = 10. FPKM values were used here as normalized expression values for correlation-based network construction and module detection, whereas differential expression inference was performed on raw counts in DESeq2. Modules were detected using a dynamic tree cut, and module eigengenes were computed. Module trait associations were evaluated with BH-adjusted p-values. Module-trait correlations were calculated using phenotypic measurements obtained from the same individual plants used for RNA-seq (n = 16). Transcription-factor enrichment was assessed using the iGRN approach, with a hypergeometric distribution test and the BH correction, to identify TFs associated with stress tolerance.

Protein interaction network analysis

Because pitaya has limited direct coverage in STRING, protein–protein interaction (PPI) networks were inferred using an orthology-based approach. Predicted pitaya protein sequences for genes in the target WGCNA modules were mapped to Arabidopsis thaliana proteins using BLASTP, retaining the best hit per gene. The interaction table was exported from the STRING protein database and visualized in Cytoscape v3.9.0^[26]. Hub nodes were ranked by degree (top 10 per module) and reported as putative conserved interaction hubs.

qRT-PCR validation

cDNA was prepared from the intended amount of RNA for each sample, based on its concentration, following all recommended protocols using the Qiagen QuantiTect Reverse Transcription Kit (Shanghai, China). Primers were designed using the pitaya coding sequences (CDS). qRT-PCR was performed using three biological replicates per treatment, with technical triplicates of each biological replicate, with the pitaya ACTIN gene (HU07G00802.1) acting as an internal control^[10]. The 2^−ΔΔCᴛ method^[27] was used to analyze the qRT-PCR results.

Statistical analysis

Statistical analysis was carried out in R. For morphological traits, Ni content, enzyme activities, and chlorophyll, values are reported as mean ± standard deviation (SD) (n = 3 biological replicates per treatment) and analyzed by ANOVA followed by Tukey's HSD at p < 0.05. For RNA-seq, two biological replicates per treatment were analyzed in DESeq2 using raw counts with BH-FDR correction; DEGs were defined at FDR < 0.05 and |log₂FC| ≥ 3.

Seedling morphology differed significantly across treatments (Fig. 1a−c). Under control (CK), seedlings showed typical growth (RL = 3.73 cm, RW = 0.12 g, SL = 8.5 cm, SW = 0.51 g, CL = 6.9 cm, CD = 5.42 mm). Melatonin alone increased biomass and modestly improved elongation traits relative to the CK. In contrast, drought (D) and nickel (Ni) significantly reduced plant height (33.49% and 40.01%), plant weight (50.05% and 59.05%), and cladode diameter (40.95% and 45.63%), respectively. Melatonin co-treatment partially restored growth under D and Ni; in DM and NiM, plant weight increased by 39.83% and 46.42% relative to D and Ni, and cladode diameter increased by 42.70% and 44.90%. Combined DNi significantly inhibited plant growth, reducing plant height, plant weight, and cladode diameter by 54.74%, 68.80%, and 53.87%, whereas melatonin supplementation (DNiM) increased these traits by 46.38%, 68.33%, and 67.46% relative to DNi (Fig. 1c). Growth traits were strongly intercorrelated (r = 0.81−0.99, p < 0.05; Fig. 1d). Melatonin altered seedling Ni accumulation in a context-dependent manner, decreasing Ni under Ni alone but increasing Ni under DNi, while still improving growth, suggesting that melatonin-mediated tolerance is not explained by reduced total seedling Ni accumulation alone and may involve detoxification and stress buffering mechanisms (Fig. 1e).

Figure 1. 

Morphological and growth responses of pitaya seedlings after 60 d under control (CK), melatonin (M), drought (D), drought + melatonin (DM), nickel (Ni), nickel + melatonin (NiM), drought + nickel (DNi), and drought + nickel + melatonin (DNiM) treatments. (a) Whole plant phenotypes. (b) Representative roots and shoots comparison. (c) Growth traits: root length (RL), root weight (RW), shoot length (SL), shoot weight (SW), cladode length (CL), and cladode diameter (CD). (d) Pearson correlation matrix among traits (*** p < 0.001; ** p < 0.01; * p < 0.05). (e) Seedling Ni accumulation (µg g^–1 dry weight). Data are means ± SD (n = 3); length in cm, weight in g, and diameter in mm. Different lowercase letters above the bars indicate significant differences among treatments at p < 0.05

Antioxidant enzymes and chlorophyll

Activities of antioxidant enzymes (APX, CAT, and SOD) increased under drought and/or Ni stress, indicating enhanced ROS-scavenging capacity (Fig. 2a−e). Melatonin generally further elevated enzyme activities in DM and NiM compared with D and Ni alone, consistent with improved tolerance. Under combined stress, DNiM showed improved growth despite APX and CAT activities slightly lower than DNi, suggesting that melatonin may rebalance redox control across multiple components rather than maximizing every single enzyme. Chlorophyll a/b declined under stress, most strongly under Ni and DNi, and melatonin partially restored chlorophyll levels under DM, NiM, and DNiM (Fig. 2a−e).

Figure 2. 

Antioxidant enzyme activities and chlorophyll contents in pitaya cladode under control (CK), melatonin (M), drought (D), drought + melatonin (DM), nickel (Ni), nickel + melatonin (NiM), drought + nickel (DNi), and drought + nickel + melatonin (DNiM) treatments. (a) Activities of ascorbate peroxidase (APX). (b) Activities of catalase (CAT). (c) Activities of superoxide dismutase (SOD). (d) Chlorophyll a content. (e) Chlorophyll b content. Data are means ± SD (n = 3). Different lowercase letters above the bars indicate significant differences among treatments at p < 0.05

Transcriptome sequencing and gene expression

We generated 16 RNA-seq libraries from pitaya cladode tissue and aligned reads to the reference genome as a quality benchmark (Supplementary Table S2). In total, 757.43 million clean reads were obtained, with mapping rates of 87.08%−95.32% (Supplementary Table S3). For exploratory quality control (QC), we computed sample-to-sample Pearson correlations and performed hierarchical clustering on FPKM values for visualization, which indicated strong concordance across replicates (r ≈ 0.859−1.0; Supplementary Fig. S1).

Differential expression across treatments

In Ni, 188 genes were differentially expressed (70 upregulated and 118 downregulated). With the application of exogenous melatonin (NiM), fewer DEGs were detected (119 total, 57 upregulated and 62 downregulated), consistent with partial normalization of Ni-responsive transcription toward the control baseline. Under D, 196 DEGs were identified (23 upregulated and 173 downregulated), whereas DM showed a shift in DEG composition (195 total, 49 upregulated and 146 downregulated). Under combined DNi stress, 167 DEGs were identified (88 downregulated and 79 upregulated), while DNiM yielded 210 DEGs (Fig. 3a; Supplementary Table S4). To specifically quantify melatonin-associated transcriptional changes under each stress background, we additionally performed within-stress contrasts, which were used for downstream pathway enrichment analysis (Pathway enrichment of DEGs section). Three DEGs (HU01G00987, HU01G00988, and novel 643) were consistently identified across all treatments, excluding the melatonin treatment (Supplementary Table S5). Analysis of log₂ fold-change (log₂FC) values showed a consistent pattern of regulation across stress conditions: HU01G00987 and HU01G00988 were consistently upregulated, while novel 643 was consistently downregulated (Fig. 3b; Supplementary Table S6). We identified TFs from 21 gene families across all stress treatments, with the AP2/ERF (Apetala2/Ethylene Response Factor), bHLH (basic Helix-Loop-Helix), and NAC (NAM, ATAF1/2, and CUC2) families being the most prominent. Most AP2/ERF TFs were downregulated under DNiM, DNi, NiM, and Ni stress, with a few upregulated under DM and DNi stress. LOB (Lateral Organ Boundaries) TFs were downregulated in DM, DNi, and D treatments. HB (Homeobox) TFs were upregulated under DNiM, DNi, and Ni stress. In contrast, SBP (SQUAMOSA Promoter Binding Protein) TFs were downregulated in both DM and D treatments. C2C2 (Zinc Finger, C2H2) TFs were upregulated in DM but downregulated in DNiM. NAC TFs were downregulated in DM but upregulated in DNi and NiM (Fig. 3c; Supplementary Table S7). We observed that AP2/ERF had the highest number of TFs under Ni and DNi stress conditions, was primarily downregulated under Ni and DNiM, and was upregulated under DM. Under DNiM, the HMG (High Mobility Group) family was most abundant, with all members downregulated. For D stress, the bHLH family had the highest number of TFs, all downregulated, while the HB was most represented under NiM stress. The bHLH TFs were mostly upregulated under Ni, NiM, and DNi. Additionally, NAC genes were significantly downregulated by more than tenfold under DM compared to CK but upregulated under DNi. Specifically, HU07G00586 from the SBP family was downregulated under both DM and D, HU05G01674 from the LOB family was downregulated under DNi, D, DM, and Ni, and HU02G00471 from the ARID family was upregulated under DM and Ni (Supplementary Table S7). The AP2/ERF, bHLH, and NAC families were the top three represented families in terms of DEGs. A heatmap of gene expression levels from these families across different treatments is provided in Supplementary Table S8 (Fig. 3d).

Figure 3. 

Differentially expressed genes (DEGs) and transcription factor (TF) families were identified under the indicated treatments. (a) Number of upregulated and downregulated DEGs under each stress treatment. (b) Heatmap illustrating the expression patterns of three representative DEGs (HU01G00987, HU01G00988, and novel 643) across treatment groups. (c) Distribution of 21 TF families represented among the identified DEGs across treatments. (d) Heatmap of the relative expression levels of 11 representative genes from the three most abundant TF families (AP2/ERF, bHLH, and NAC) across treatment groups. DEGs were defined at FDR < 0.05 and |log₂FC| ≥ 3.

Gene ontology enrichment analysis

GO enrichment analysis showed that DEGs with distinct expression patterns were enriched across various categories (Supplementary Table S9). Upregulated DEGs exhibited consistent enrichments under single and combined stress treatments; therefore, we analyzed the top five GO pathways. Generation of precursor metabolites and energy (GO:0006091), drug metabolic process (GO:0017144), carboxylic acid metabolic process (GO:0019752), organic acid metabolic process (GO:0006082), and small molecule metabolic process (GO:0044281) were consistently enriched across all treatments. Genes HU01G00987 and HU01G00988 co-accumulated and participated in these pathways, and HU03G02702 was involved across all stress conditions except D stress (Fig. 4; Supplementary Table S10).

Figure 4. 

Gene Ontology (GO) enrichment analysis of DEGs under different treatments. Panels (a)−(g) correspond to a specific stress treatment: (a) melatonin; (b) drought; (c) drought + melatonin; (d) nickel; (e) nickel + melatonin; (f) drought + nickel; (g) drought + nickel + melatonin. For each treatment, the x-axis shows the gene ratio, and the y-axis lists enriched GO terms. Dot size indicates the number of DEGs, and dot color scales with FDR (darker colors indicate more substantial enrichment).

For DNi and NiM, the top five GO pathways involving downregulated genes were in the MF category. Two binding-related pathways, heme binding (GO:0020037) and tetrapyrrole binding (GO:0046906), showed consistent co-enrichment across Ni, DNi, NiM, and DM. The pathways with the highest gene enrichment across all stress treatments were hydrolase activity acting on glycosyl bonds (GO:0016798) and hydrolase activity hydrolyzing O-glycosyl compounds (GO:0004553). Notably, HU04G00837, which functions in both hydrolase pathways, was enriched under DNi, DNiM, and Ni, but not under D or NiM, consistent with the pattern for upregulated genes in which D stress showed distinct pathway enrichment profiles (Supplementary Table S11). This indicates that the pathways activated under D stress are significantly different from those triggered by the other treatments.

Under DNiM treatment, non-membrane-bounded organelles and intracellular non-membrane-bounded organelles were significantly enriched. Hydrolase activity acting on a glycosidic bond and protein heterodimerization activity also showed enrichment across DNiM. Under DNi stress, the generation of precursor metabolites and energy, phosphoenolpyruvate carboxykinase activity, and phosphoenolpyruvate carboxylase activity were significantly enriched. Hydrolase activities acting on glycosyl bonds, carbon-carbon lyase, and carboxy-lyase activities, and N-glycosyl compounds likewise showed significant enrichment.

Under NiM stress, pathways related to energy derivation by oxidation of organic compounds, tricarboxylic acid metabolism/cycle, aerobic respiration, and carbon fixation were enhanced. Similarly, Ni stress-enriched pathways are associated with aerobic respiration, the tricarboxylic acid cycle, and carbon fixation. Under D stress, hydrolase activity, hydrolysis of compounds, and protein dimerization activity were enriched. Consistently, the functional categories present in each treatment exhibited partially overlapping enrichment patterns.

To connect enrichment results to candidate genes, we examined DEGs within representative significant GO pathways (p-adj < 0.05) across treatments (Supplementary Table S12). For example, hydrolase activity, acting on glycosyl bonds (GO:0016798), contained 22 genes that were predominantly downregulated, with a subset showing treatment-dependent patterns. Likewise, protein heterodimerization and N-glycosyl-compound hydrolase activities (GO:0046982, GO:0016799) were enriched mainly among downregulated DEGs, particularly under D and DNiM. In contrast, carbon fixation and phosphoenolpyruvate carboxykinase activity (GO:0015977, GO:0004611, GO:0008964) shared common upregulated genes across multiple stress conditions (Supplementary Table S12).

Pathway enrichment of DEGs

To characterize pathways enriched among DEGs in within-stress contrasts that isolate melatonin application under D and Ni exposure, we conducted KEGG pathway enrichment analysis.

In DM vs D, we found 20 significantly enriched pathways. DEGs associated with melatonin application were enriched in carbon/specialized metabolism and stress signaling pathways. Enriched terms included phenylpropanoid and flavonoid biosynthesis, multiple alkaloid routes (indole, isoquinoline, tropane/piperidine/pyridine), amino sugar and nucleotide sugar metabolism, phenylalanine/tyrosine/tryptophan biosynthesis, alpha-linolenic acid metabolism, and biosynthesis of unsaturated fatty acids. Stress signaling categories, such as plant MAPK signaling and circadian rhythms, were significant in plants. At the transcriptome level, enrichment patterns were consistent with increased allocation of carbon toward energy/redox support. Collectively, these enrichment patterns were consistent with diversion of carbon into antioxidant/structural phenolics and lipid-derived oxylipins, alongside shifts in signaling pathways (e.g., MAPK signaling and circadian rhythm terms) associated with drought responses.

In NiM vs Ni, we found nine significantly enriched pathways. DEGs associated with melatonin application were enriched in core energy and detox pathways, pyruvate metabolism, fatty acid metabolism/biosynthesis, and propanoate metabolism, as well as glutathione metabolism, which support ROS quenching and metal chelation. We also observed enrichment in porphyrin metabolism, consistent with partial restoration of chlorophyll/photochemistry, and translational capacity terms (ribosomes, spliceosome), suggesting enhanced protein synthesis for repair and defense.

In DNiM vs DNi, we found 16 significantly enriched pathways. With both stresses imposed, melatonin-associated DEGs additionally showed enrichment for membrane trafficking and barrier formation pathways, SNARE interactions in vesicular transport, phagosome formation, motor proteins, and cutin/suberin/wax biosynthesis. Central metabolism and redox hubs are also significant, including carbon metabolism, the pentose phosphate pathway, and peroxisome and cysteine/methionine metabolism. Plant hormone signal transduction, MAPK signaling, and phenylpropanoid biosynthesis completed the response, indicating coordinated reinforcement of membranes, trafficking, and antioxidant capacity under compound stress.

Regulation of ROS and cell wall-related genes

Our pathway analysis revealed differential regulation: genes within MAPK signaling and carotenoid biosynthesis pathways were downregulated, while those in the pentose phosphate pathway and carbon metabolism were upregulated. Notably, genes active in carbon metabolism under Ni exposure persisted under DM conditions, with even lower transcript abundance. HU11G01285, a gene annotated to phenylpropanoid biosynthesis, exhibited consistent expression across Ni, DNi, and NiM treatments.

Under NiM, genes involved in the pentose phosphate pathway and glutathione metabolism were transcriptionally upregulated. In contrast, genes involved in carotenoid and phenylpropanoid biosynthesis were downregulated, and ultimately, genes linked to polysaccharide and cell wall metabolism were also downregulated. In this study, we found that most of the related genes were enriched in amino sugar and nucleotide sugar metabolism pathways, and all were uniformly downregulated. Notably, specific carbon metabolism genes (HU01G00987, HU03G02702, HU01G00988, HU10G01890) showed transcript expression under both NiM and DM conditions, with reduced transcript abundance in DM treatment (Supplementary Table S13).

Transporters and metabolic pathways

Transport-associated genes were transcriptionally upregulated specifically under DM and M. The DM treatment showed the highest functional enrichment for genes related to cellular transport, including HU01G00987, HU03G02702, HU01G00988, HU01G00991, and HU10G01890, all of which were annotated to carbon metabolism pathways. Additionally, HU08G01993 was upregulated and annotated to starch and sucrose metabolism, while HU05G00487 and HU08G01254 were transcriptionally downregulated.

Energy metabolism pathways

Our results revealed that DM treatment showed significantly enriched DEGs related to energy metabolism, with all genes involved in carbon metabolism, pyruvate metabolism, and carbon fixation in photosynthesis showing transcriptional upregulation. Similarly, D treatment triggered a significant number of DEGs related to energy metabolism, characterized by transcriptional upregulation of all pyruvate metabolism genes. In contrast, genes linked to glycerolipid metabolism were transcriptionally downregulated.

Transcription factors and stress pathways

In the DM treatment, stress-responsive DEGs showed consistent downregulation across multiple pathways, including cutin, suberin, and wax biosynthesis; plant MAPK signaling; and glutathione metabolism. Similarly, D treatment induced downregulation of genes (HU01G00570, HU07G01997, HU03G00839, HU07G00657, HU10G00264, HU02G03030, HU09G01540, HU06G01980, HU09G00700, HU03G00194, HU01G00602) involved in plant hormone signal transduction, carotenoid biosynthesis, phenylpropanoid biosynthesis, cutin, suberin, and wax biosynthesis, and plant MAPK signaling pathways. Notably, genes involved in carotenoid biosynthesis were transcriptionally downregulated in the Ni (HU02G01579), DNi (HU02G01579, HU04G00678), NiM (HU02G01579), DNiM (HU02G01579), and D (HU02G01579) treatment groups.

Photosynthesis and carbohydrate metabolism

Our results revealed that DNi treatment explicitly regulated genes involved in carbon fixation in photosynthetic organisms and in the pentose phosphate pathway. In addition, genes involved in amino sugar and nucleotide sugar metabolism were transcriptionally downregulated. The NiM treatment consistently transcriptionally upregulated genes associated with carbon metabolism, the pentose phosphate pathway, and photosynthetic carbon fixation pathways.

WGCNA modules and trait associations

Network construction

Across the 16 RNA-seq libraries, 30,793 genes were expressed after filtering low-abundance transcripts. Using pickSoftThreshold, β = 10 yielded a scale-free topology fit (R² > 0.8), with acceptable mean connectivity, and was used for the unsigned network. Accordingly, the dynamic tree cut with merging yielded 21 gene co-expression modules (Fig. 5), with module sizes ranging from 26 to 273 genes. Module eigengenes (MEs) summarized each module's expression profile. Genes in several modules showed coherent regulation across treatments. MEblack, MEblue, and MEgreen were predominantly transcriptionally downregulated, while MEgreenyellow, MElightcyan, MElightgreen, MEroyalblue, and MEsalmon were predominantly transcriptionally upregulated (Supplementary Tables S12, S14). Module trait correlations (RL, RW, SL, SW, CL, CD) were BH-adjusted, and FDR-significant associations are summarized in the supplementary material.

Figure 5. 

WGCNA-based co-expression modules and transcription factor (TF) distribution. (a) Effect of soft-thresholding power on scale-free topology fit and mean connectivity. The red horizontal line indicates the scale-free topology fit threshold of R² = 0.80 used to guide soft-thresholding power selection for WGCNA network construction. (b) Hierarchical clustering dendrogram of expressed genes based on topological overlap, showing assignment of genes into co-expression modules identified by WGCNA. TFs were mainly concentrated in C2C2, MYB, C2H2, HB, C3H, bHLH, AP2/ERF, bZIP, and WRKY families, consistent with their potential roles in abiotic stress regulation.

Identification of stress-associated modules

Among the 21 gene modules, three exhibited higher specificity for stress tolerance. Based on the heatmap of the module-trait correlations (Fig. 6b), the MEgreenyellow module showed the highest negative correlation with RL (r = –0.842, p = 4.2 × 10^–5). In contrast, the MEgreen module exhibited the highest positive correlation (r = 0.729, p = 0.0014). In addition, the MEpink module had the highest correlation with SL (r = 0.615, p = 0.011) (Fig. 6a). The MEgreen and MEpink modules showed a significant positive correlation, while the MEgreenyellow module showed a negative correlation. This indicates that the genes in the MEgreen and MEpink modules may share similar functions or regulatory roles under certain conditions. In contrast, the MEgreenyellow module may be linked to opposing biological processes. The gene expression profiles across the three modules showed significant phenotype-specific variation, suggesting that specific genes may play a crucial role in regulating a particular trait. We screened these modules for further analysis to identify candidate genes that play significant roles in RL and SL of pitaya under single and combined stresses. By identifying protein–protein interactions (PPIs), we selected the top 10-degree hub genes from the network (Supplementary Table S15). The network interactions of hub genes within the three modules are shown in Fig. 6c. Because the STRING interactions were retrieved from Arabidopsis thaliana based on orthology, the networks in Fig. 6c represent putative conserved PPIs rather than experimentally validated pitaya-specific interactions.

Figure 6. 

Module-trait relationship and hub-gene networks under abiotic stress. (a) Scatterplots showing the relationships between module membership and gene significance for the three selected modules (MEgreen, MEpink, and MEgreenyellow). Module membership indicates the correlation between the expression of an individual gene and the corresponding module eigengene, whereas gene significance indicates the correlation between that gene and the associated phenotypic trait. The x-axis represents module membership, while the y-axis displays gene significance for traits. (b) Full heatmap of correlation between 21 WGCNA modules and six phenotypic traits: shoot length (SL), cladode length (CL), root length (RL), root weight (RW), shoot weight (SW), and cladode diameter (CD). Red indicates positive correlation and blue indicates negative correlation. Each cell reports the correlation coefficient and the BH-adjusted p-value (FDR). (c) Protein–protein interaction (PPI) network of the top 10 hub genes identified in the MEpink, MEgreen, and MEgreenyellow modules. Nodes represent hub genes, and edges indicate predicted conserved interactions. These networks were inferred from Arabidopsis orthologs using STRING and therefore represent orthology-based predicted interactions rather than experimentally validated pitaya-specific PPIs.

Representative hub genes within the MEpink, MEgreen, and MEgreenyellow modules were associated with pathways including ubiquitin-mediated proteolysis, amino sugar and nucleotide sugar metabolism, diterpenoid biosynthesis, and carotenoid biosynthesis. They are also associated with Gene Ontology (GO) biological processes, including regulation of protein metabolic process, cell cycle, regulation of cellular protein metabolic process, monosaccharide metabolic process, regulation of molecular function, regulation of catalytic activity, regulation of phosphorus metabolic process, regulation of protein modification process, and regulation of phosphate metabolic process.

In the GO analysis, the MEgreen module was significantly enriched for nucleus and organelle parts, intracellular organelle parts, and the spliceosome. Biosynthesis of amino acids and protein processing in the endoplasmic reticulum were highlighted in the KEGG analysis. Meanwhile, the MEpink module in NiM was enriched for transcription regulator activity and DNA-binding transcription factor activity. It showed significance in plant–pathogen interaction and plant hormone signal transduction in KEGG analyses. In analyses of heme binding, tetrapyrrole binding, and plant hormone signal transduction, the MEgreenyellow module was highlighted. Among them, the MEgreen module was the main target of the terms related to 'organelle part' in the GO analysis, suggesting that this module is closely linked to the specific functions of organelles, such as protein synthesis, energy production, and material transport. Additionally, the strong enrichment of plant-pathogen interaction and plant hormone signal transduction pathways within these modules highlights their key role in regulating plant growth, development, and responses to external environmental stimuli (Fig. 7; Supplementary Tables S16, S17).

Figure 7. 

Relative expression patterns of genes in representative co-expression modules. Heatmaps show the relative expression profiles of genes assigned to the three representative modules: (a) MEgreen; (b) MEpink; (c) MEgreenyellow module. Rows represent genes within each module, and the columns represent treatment samples. The x-axis shows the treatment groups in order: control (CK); melatonin (M); drought (D); drought + melatonin (DM); nickel (Ni); nickel + melatonin (NiM); drought + nickel (DNi); drought + nickel + melatonin (DNiM); with biological replicates one and two for each treatment (n = 2). Color intensity represents relative normalized expression across samples, illustrating treatment-associated expression variation within each module.

qRT-PCR validation

To validate the gene expression results obtained from RNA-seq, eight candidate genes (HU01G00987, HU01G00988, HU04G00837, HU01G00991, HU08G02052, HU04G00898, HU04G00678, HU11G01840) under various treatments were selected for qRT-PCR analysis (Supplementary Table S18). We compared the 2^−ΔΔCᴛ values from qRT-PCR (SD < 0.5) with the log₂FC values from RNA-seq for each treatment group. The results showed that the expression profiles of most candidate genes across different treatments were consistent with the RNA-seq data, confirming the reliability of the analyses conducted in this study (Fig. 8). The validated genes collectively represent biologically relevant stress-response functions, including central carbon metabolism (HU01G00987, HU01G00988, HU01G00991), cell wall remodeling (HU04G00837), and hormone/stress signaling (HU04G00678), indicating that the qRT-PCR-supported candidates span multiple processes associated with pitaya responses to drought, nickel, and combined stress.

Figure 8. 

qRT-PCR validation of eight candidate genes and comparison with RNA-seq expression patterns. Blue bars represent relative qRT-PCR expression levels, and orange lines represent the RNA-seq expression values across CK, M, D, DM, Ni, NiM, DNi, and DNiM treatments. Error bars indicate the standard deviation of RT-qPCR measurements from three biological replicates. The selected genes were HU01G00987.1, HU01G00988.1, HU04G00837.1, HU01G00991.1, HU08G02052.1, HU04G00898.1, HU04G00678.1, and HU11G01840.1.

Climate change and human activities are driving more frequent drought and heavy metal contamination, which together threaten crop productivity and global food security by disrupting physiological, biochemical, and molecular processes. Heavy metals such as nickel can accumulate in agricultural soils, constrain crop productivity, and complicate the safe utilization of harvested biomass and downstream processing streams^[28]. Pitaya, a member of the cactus family, is generally considered tolerant to drought and heat^[29], but the pitaya responses to single and combined D and Ni stress, and how melatonin may mitigate these stress conditions, remain poorly understood. Further elucidation of these responses may help develop pitaya varieties that are resilient to both D and heavy metal stress. Additionally, elevated Ni levels in pitaya can negatively affect the quality and marketability^[30]. In the context of sustainable agriculture and fruit peel recycling, excessive Ni content could compromise product safety and commercial value^[31]. Prior studies report that melatonin can reduce or increase heavy metal accumulation, depending on species, tissue, and stress context^[32,33]. In our research, melatonin's effect on seedling Ni content was context-dependent, consistent with tolerance mechanisms rather than a uniform reduction in seedling Ni accumulation. A similar pattern has been reported in pitaya under combined heavy-metal and environmental stress, where melatonin-treated plants accumulated higher shoot metal content yet still showed improved growth, supporting the view that stress mitigation may involve tissue-level physiological tolerance and detoxification rather than simple metal exclusion alone^[2,10,21,34].

In our study, exogenous melatonin application significantly promoted pitaya seedling growth across all stress treatments, particularly under combined DNi stress, where melatonin partially restored height, biomass, and cladode diameter compared with DNi alone. These findings are consistent with previous reports that melatonin enhances growth under stress conditions^[35]. Based on phenotypic data, pitaya seedlings appeared more sensitive to Ni than to D, suggesting that pitaya plants may have higher tolerance to D than to Ni. Because pitaya exhibits CAM-like physiology, the PEG treatment used here should be interpreted as a controlled osmotic drought-like treatment rather than a complete representation of field drought conditions in a cactus crop; accordingly, our drought-related conclusions are centered on responses under the imposed experimental system.

Antioxidant enzymes and chlorophyll under stress

Antioxidant enzymes play a crucial role in plant defense against oxidative stress. Stress-induced ROS accumulation is a common feature of plant responses to abiotic stress and contributes to redox imbalance under adverse conditions^[36,37]. In our study, activities of APX, CAT, and SOD increased under drought and nickel stress, indicating enhanced ROS-scavenging capacity, particularly for H₂O₂, which aligns with previous findings^[38,39]. In contrast, chlorophyll a and b declined, especially under Ni and DNi, consistent with Ni interfering with chlorophyll biosynthesis and nutrient uptake, such as iron and manganese. The more substantial reduction of chlorophyll under DNi than under either single stress is consistent with a greater physiological burden under the combined stress treatment. Generally, melatonin is reported to boost antioxidant capacity and protect the photosynthetic machinery^[40]. In our experiment, melatonin increased APX, CAT, and SOD under most stress conditions and partially restored chlorophyll levels. However, under DNiM, APX and CAT activities were slightly lower than under DNi, despite improved growth. This pattern may reflect redistribution of redox regulation across multiple enzymatic and non-enzymatic components, and suggests that improved tolerance does not necessarily require maximal activity of each antioxidant enzyme. This physiological pattern is also consistent with the transcriptomic data. Melatonin-associated treatments showed enrichment of redox-related pathways, including glutathione metabolism and peroxisome-associated functions, as well as porphyrin metabolism under NiM. These changes provide a transcript-level context for the observed antioxidant responses and partial chlorophyll recovery.

Gene expression triggered by D and Ni stress

Under drought stress, plant cell walls can be flexibly remodeled. Additionally, nickel has been shown to reduce cell wall plasticity and alter cell wall chemical composition^[41]. The toxicity of Ni in plant stems disrupts the ionic balance and damages cell membranes. In our study, we observed the upregulation of phosphoenolpyruvate carboxylase (PEPC), α-amylase, and the transcription factor bHLH147 under drought stress. Notably, bHLH TFs have been shown to play a crucial role in regulating the synthesis of cell wall components, including cellulose, hemicellulose, and pectin^[42]. Under Ni stress, proteins highly associated with cell wall synthesis, including expansin-like B1, α-glucan water dikinase 1 (GWD1), homeobox leucine zipper protein ATHB-12, bidirectional sugar transporter SWEET17, and scopoletin glucosyltransferase, showed upregulation. This expression pattern may reflect Ni-associated stress responses affecting cell wall-related processes. Additionally, under Ni stress, the gene HU06G02339, which encodes CSLG1, and the gene HU01G00606, which encodes the leaf rust 10 disease-resistance locus receptor-like protein kinase 1.2 (LRK10L1.2), were upregulated.

Under D treatment, we observed upregulation of the gene HU01G02786, which encodes pyruvate kinase isozyme A, and of the gene HU02G03020, which encodes omega-amidase. Additionally, HU08G00833 (sp|F4JJR8|FEY_ARATH), HU05G01295 (sp|Q9LHI0|NDUA6), HU01G00988 (phosphoenolpyruvate carboxylase), HU0G00205 (α-glucan water dikinase 1, GWD1), HU01G00157 (scopoletin glucosyltransferase), and HU11G01114 (alcohol dehydrogenase) were also upregulated. Together, these expression patterns are consistent with adjustments in cell wall-related processes and central carbon metabolism under Ni and D stress. These transcript-level changes also provide a plausible molecular context for the observed reductions in growth and cladode diameter under Ni and D treatments. They may also help explain the partial recovery of biomass-related traits following melatonin application.

TFs identified under D and Ni stresses

Several transcription factors were highlighted, including the homeobox-leucine zipper protein ATHB-12, MYB102, bHLH147, and bHLH30, as well as AtHsp70-1 and a Zinc finger HIT domain-containing protein, suggesting engagement of conserved regulatory nodes under D and Ni stress in pitaya^[43]. MYB transcription factors, particularly the R2R3-MYB group, are widely implicated in abiotic stress responses by modulating secondary metabolism, ROS homeostasis, and hormone signaling^[44−47]. Likewise, bHLH transcription factors contribute to drought and heavy-metal stress tolerance by regulating stomatal behavior, photosynthesis, and nutrient/metal homeostasis^[47−53]. The presence of MYB and bHLH candidates in our dataset is therefore consistent with possible TF-mediated integration of metal and drought signals, linking upstream stress to inputs to downstream metabolic and redox responses. Such TF-associated regulation may help explain the better growth and physiological performance of melatonin-treated plants under stress, despite substantial transcriptomic perturbation.

Impact of stress treatments on key metabolic pathways and gene expression

At the transcriptome level, we retained 30,793 expressed genes in pitaya seedlings, and most stress-responsive DEGs are enriched in pathways related to polysaccharide and cell wall metabolism, ROS detoxification, cellular transport, signal transduction, energy metabolism, responses to abiotic stress, biosynthesis of secondary metabolites, carbohydrate metabolism, and antioxidant systems^[54]. These processes are central to maintaining cellular homeostasis and structural integrity, regulating metabolic networks, and enhancing metabolic flexibility and adaptability under stress^[55]. Such metabolic adjustments are closely linked to ATP production, redox homeostasis, and biosynthetic precursor supply in plant tissues under stress^[56]. Genes involved in carbon metabolism were broadly upregulated in M, DM, Ni, NiM, and DNi treatments, whereas drought alone produced a mixed pattern, with some genes upregulated and others downregulated. For genes associated with starch and sucrose metabolism, expression was downregulated in both melatonin and Ni stress treatments, consistent with reports that stress can reduce cantaloupe sweetness through altered sucrose accumulation. Together with the chlorophyll and growth data, these results are consistent with a melatonin-associated shift in carbon-energy-redox pathways that may help sustain growth under stress.

To improve the biological interpretation of the trait-linked modules, we further examined representative hub genes using their putative ortholog annotations. In the MEgreenyellow module, HU11G01139 showed similarity to Arabidopsis SKP2B, an F-box protein associated with SCF-type ubiquitin ligase complexes^[57], suggesting that stress-responsive protein turnover may contribute to growth adjustment under severe stress. In the same module, HU04G01748 was annotated as ATHB-12-like, a homeobox-leucine zipper transcription factor linked to abiotic stress responses^[58], indicating a possible role in coordinating developmental and stress signaling responses. In the MEpink module, HU04G00678 showed similarity to CYP707A2, an ABA 8-hydroxylase, suggesting potential involvement in ABA turnover and drought-associated signaling^[59]. In the MEgreen module, HU08G00128 was annotated as an early nodulin-like protein, while HU04G00837 was annotated as alpha-L-arabinofuranosidase 2-like, suggesting possible roles in extracellular signaling and cell wall remodeling that may influence elongation-related traits. In addition, HU09G01329 encodes a BURP-domain protein-like sequence, which may be associated with stress-responsive structural adjustment^[60]. Together, these putative functions suggest that the hub genes may be associated with protein turnover, hormone regulation, and structural remodeling during pitaya acclimation to drought, nickel, and combined stress, although these assignments remain provisional and require functional validation.

This study focused on morphophysiological and transcriptomic responses to exogenous melatonin under drought, nickel, and combined stress in pitaya. This study did not include target validation, receptor antagonists/structural analogues, or genetic perturbations; therefore, our transcriptomic and network results should be interpreted as associations rather than evidence of receptor-dependent melatonin signaling. Because RNA-seq was conducted with two biological replicates per treatment, transcriptomic patterns were interpreted together with stringent DEG filtering, replicate concordance, independent morphophysiological observations, and qRT-PCR validation of representative genes. The data are consistent with melatonin-associated modulation of carbon-energy metabolism, redox buffering, and growth-related regulatory programs under D, Ni, and DNi stress.

Exogenous melatonin enhanced pitaya performance under drought, nickel, and especially the combined drought-nickel stress conditions, mainly by coordinating carbon-energy-redox adjustments rather than by significant reductions in tissue nickel accumulation. Under D and Ni stress, growth declined relative to the control, whereas melatonin partially improved performance together with antioxidant activity and chlorophyll retention. Transcriptome profiling indicated melatonin-associated enrichment of central carbon flow (carbon metabolism, starch/sucrose turnover, and pyruvate-linked routes), redox buffering (peroxisome/glutathione), and signaling layers (MAPK and hormone pathways). Co-expression analysis resolved 21 modules; MEgreen and MEpink were positively associated, whereas MEgreenyellow was negatively associated with elongation traits. Protein–protein interaction (PPI) mapping within trait-linked modules highlighted hub genes associated with ubiquitin-mediated proteolysis and carbohydrate/secondary metabolism, and qRT-PCR of representative candidates corroborated RNA-seq trends. Collectively, these results are consistent with melatonin-associated modulation of metabolic supply and redox buffering under imposed stress conditions. The most biologically relevant candidate genes HU11G01139 (SKP2B), HU04G01748 (ATHB-12-like), HU04G00678 (CYP707A2-like), HU08G00128 (ENODL1), HU04G00837 (alpha-L-arabinofuranosidase 2-like), and HU09G01329 provide targets for functional validation, bio-stimulant optimization, and breeding to improve pitaya resilience.