Effects of Organic Selenium and Nano-Selenium on Drought Stress of Pak Choi (Brassica chinensis var. pekinensis. cv. ‘Suzhouqing’) and Its Transcriptomic Analysis

1. Introduction

‘Suzhouqing’ (Brassica chinensis var. pekinensis) plays a crucial role in the annual vegetable supply due to its low crude fiber content, delightful texture, and abundant nutritional value [1]. However, the frequent occurrence of severe drought in China has emerged as a significant constraint on the growth of Suzhouqing. Drought stress directly weakens photosynthetic intensity, and triggers excessive production of reactive oxygen species within cells, disrupting the balance of oxidative metabolism and leading to intensified membrane lipid peroxidation damage or even plant mortality [2,3,4].

Selenium is an important element for crop growth and plays a vital regulatory role in crop development. It actively participates in the formation of antioxidant enzymes (GSH-PX) and facilitates the conversion of harmful peroxides into harmless hydroxyl compounds, thereby safeguarding biofilms against oxidative degradation [5]. Exogenous application of selenium is a key approach to enhance crop drought resistance, which not only increases the selenium content in the entire food chain but also improves crop yield and quality. Previous studies have demonstrated that selenium can promote crop growth rate, mitigate oxidative damage caused by ultraviolet radiation, increase chlorophyll and carotenoid content in plant leaves, as well as boost antioxidant enzyme activity and osmoregulatory substance content [6,7,8,9]. It also contributes to increased selenium content and yield of fruits, effectively enhancing fruit quality and improving stress resistance in crops [10]. Moreover, research has shown that applying selenium fertilizer under drought stress effectively improves plant tolerance to drought. The detrimental impacts of drought stress on wheat [11], barley [12], and rape [13] were alleviated through enhanced photosynthesis capacity, increased antioxidant capacity, and elevated osmoregulatory substance content in plants.

However, there is limited research on the regulation of exogenous selenium on the physiological, biochemical, and metabolic growth of pak choi under drought stress. Most studies focusing on the effects of selenium on plant growth and nutritional quality primarily measure physiological indicators, with few reports considering the molecular mechanisms involved. This limitation hampers our comprehensive understanding of the physiological function and mechanism of selenium in plants and impedes the utilization of selenium-rich plant resources and development of selenium-enriched agricultural products. Therefore, photosynthetic pigment, antioxidant, nutritional quality and other physiological indexes were determined, and the molecular mechanism of the effects of organic selenium and nano-selenium on leaves was analyzed by transcriptomics in this study.

2. Materials and methods

2.1. Experimental Site and Experimental design

The test material was Suzhouqing, and the seeds were purchased from Jiangxi Ruibao Seed Co., LTD (Jiujiang, China). Nano-selenium, a tiny red elemental nanoparticle of selenium with a particle diameter ranging from 20 to 60 nm, was developed by the Zhongnong Selenium-rich Agricultural Technology Research Institute in Beijing, China. On the other hand, organic selenium, mainly consisting of malt selenium ((CH3)2Se), was developed by Northwest A&F University in Yangling, China.

The experiment was conducted from November 2022 to January 2023 in the laboratory of the College of Land Resources and Environment, Jiangxi Agricultural University, Nanchang, Jiangxi Province, China. After 10 days of sowing, the seedlings were selected and transplanted into pots (33.00 × 24.40 × 13.50 cm) with 2 plants per pot. The medium was organotrophic soil, and its physical and chemical properties are shown in Table 1.

When the seedlings grow to the four-leaf stage (December 22nd, 2022), natural drought was adopted to control water, and three pots were repeated for each treatment. We used an electronic scale with 0.50 g sensitivity to weigh at 18:00 every day, and evenly supplemented the lost water on the day, controlled the soil moisture content at (10.00±0.50)%, and recorded the water consumption.

Three levels of exogenous selenium (20 mL) were applied, including nano-selenium treatment: 5, 10, 20 mg L⁻¹ (DN5, DN10, DN20), organic selenium treatment 5, 10, 20 mg L⁻¹ (DO5, DO10, DO20), and a control group (DCK: 20 mL water). Foliar spraying was conducted on December 9th and 16th, 2022 (before drought treatment), as well as on December 27th, 2022 (during drought treatment).

After 10 days of drought treatment (January 1st, 2023), fresh leaves were taken, part of them were immediately tested for relevant physiological indexes in vitro, and part of the samples were placed in the refrigerator at -80 °C.

2.2. Determination of physiological and biochemical indexes

The concentrations of chlorophyll a, b and carotenoids were determined using an ultraviolet spectrophotometer (UV-2450, Shimadzu, Japan) at wavelengths of 665 nm, 649 nm, and 470 nm respectively. Acetone-ethanol mixed extraction was employed for the analysis, and the approach was cited in the publication authored by Wang and Huang [14]. The LCpro T portable photosynthetic measurement system was utilized to determine various photosynthetic indices including net photosynthetic rate (Pn), stomatal conductance (Gs), transpiration rate (Tr), and intercellular CO₂ concentration (Ci) during the peak functional period of leaves from 9:00 am to 11:00 am on sunny days.

The activity of superoxide dismutase (SOD) was determined using a novel highly water-soluble tetrazolium salt, WST-8 (4-[3-(2-methoxy-4-nitrophenyl)-2-(4-nitrophenyl)-2h-5-tetrazolium]-1, 3-benzenesulfonate), as described by Ukeda et al. in 1999 [15]. Superoxide ions generated by xanthine oxidase were able to reduce WST-8 to water-soluble formaldehyde, resulting in maximum absorbance at 460 nm.

The activity of peroxidase (POD) was assessed by monitoring the oxidation rate of guaiacol. The reaction was initiated by adding 200 μL of enzyme solution to a centrifuge tube containing 10 mM H₂O₂, 25 mM sodium phosphate buffer and 0.05% guaiacol. After thorough mixing, the change in absorbance at 470 nm was measured over a period of 90 seconds, and subsequently used to calculate the POD activity [16].

We utilized the peroxidatic function of catalase to determine catalase (CAT) activity. The method was based on the enzymatic reaction with methanol in the presence of an optimal concentration of hydrogen peroxide [17]. The resulting formaldehyde was quantified spectrophotometrically using 4-amino-3-hydrazino-5-mercapto-1,2,4-triazole (Purpald) as a chromogen.

The malondialdehyde (MDA) content was determined using the thiobarbituric acid reaction method [18]. We measured 0.10 g of the sample, added 1 mL of the extract, and homogenized it at room temperature. Subsequently, we centrifuged the mixture at 10,000 rpm for 10 min at room temperature to obtain a supernatant volume of 100 μL. This supernatant was then combined with 300 μL of malondialdehyde (MDA) content kit reagent and incubated in tightly covered water at a temperature of 92 ℃ for 30 min to prevent water loss. After cooling down by centrifuging again at 12,000 rpm for another 10 min at a temperature of 25 ℃, we transferred an aliquot containing 200 μL supernatant into a well on a standard microplate and measured its absorbance values A₅₃₂ and A₆₀₀. The difference between these two absorbance values was denoted as ΔA = A₅₃₂-A₆₀₀.

Precisely measured by Nishimoto et al., the content of reduced glutathione (GSH) was observed [19]. Approximately 0.10 g of samples were weighed and 1 mL of extraction liquid was added. Following homogenization in an ice bath, the resulting supernatant was centrifuged at a speed of 12,000 rpm and a temperature of 4 ℃ for a duration of 10 min. The obtained supernatant was then transferred to an ice-filled container for subsequent measurement. Distilled water, superserum, and two reagents from the reduced glutathione (GSH) kit were introduced into a 96-well plate, thoroughly mixed, and allowed to stand for a period of 2 min. Absorbance values at a wavelength of 412 nm were recorded for both the assay tube and blank tube. Finally, the GSH content was determined based on the difference in absorbance values.

The determination of soluble sugar content was conducted following the method described by Buysse and Merckx [20]. Briefly, 0.10 g of leaf sample was weighed and mixed with 1 mL of 80% ethanol. The resulting supernatant was then heated in a water bath at 80 ℃ for 30 min. This process was repeated twice, and the combined supernatants were used as the extraction solution. Subsequently, a volume of 50 μL extract was diluted by a factor of 40 and mixed with 5 mL anthrone sulfuric acid reagent. After boiling for 10 min, the mixture was removed from heat and cooled to room temperature using tap water before measuring its absorbance at a wavelength of 620 nm. Finally, the total soluble sugar content was determined based on a standard curve.

The Bicinchoninic acid (BCA) assay, which involved two reactions, was employed for the determination of soluble protein content [21]. In the initial reaction, a copper ion complex with amide bonds was formed, resulting in the reduction of copper in an alkaline solution. Subsequently, the BCA reagent was primarily reduced by the copper–amide bond complex and also by tyrosine and tryptophan residues.

The amino acid content was determined according to the method described by Chen et al. [22]. Approximately 0.10 g of leaves were extracted with 10 mL of distilled water in a boiling water bath for 45 min. After centrifugation for 10 min, 1 mL of supernatant was collected and mixed with 0.5 mL of a 2% ninhydrin reagent and 0.5 mL of pH 8.0 phosphate buffer solution. The reaction mixture was heated in a boiling-water bath for 15 min and then cooled to room temperature before adding an additional volume of distilled water (8 mL). After standing for another 10 min, the absorbance at a wavelength of 562 nm was measured using a MetashUV-5200 UV-vis spectrophotometer.

2.3. RNA extraction and RNA-Seq sequencing

Plant total RNA was extracted using a specialized kit. Quality assessment involved agarose gel electrophoresis to determine degradation levels and contamination status, while RNA concentration was quantified using Qubit2.0 fluorometer (Thermo Scientific, USA) technology. The integrity of the RNA samples was precisely evaluated with the Agilent 2100 system (Agilent Technologies Co. Ltd., USA) before proceeding to library construction for subsequent Illumina sequencing, which followed Nanodrop analysis to determine sample purity [23].

2.4. Identification of gene differential expression

To ensure the quality and reliability of data analysis, we utilized Trimmomatic 0.39 software (version: Trimmomatic 0.39) to filter the sequencing data. Specifically, we removed a subset of reads from the original data that contained adapters and exhibited low sequencing quality [24]. Subsequently, Trinity was employed to concatenate the clean reads, resulting in the generation of Unigene sequences which served as reference sequences (Ref) [25]. RSEM software was then applied to compare each sample’s clean reads with Ref. Any reads with comparative quality lower than 10 were filtered out accordingly. Additionally, any unpaired reads were aligned to multiple regions of the genome. Finally, RSEM v1.2.28 software was used to calculate the TPM value (Transcripts Per Kilobase Million) for each sample in order to estimate gene expression levels. Ultimately, genes exhibiting significant differences in expression levels between nano-selenium and organic selenium were identified using statistical methods based on criteria q-value < 0.05 [23].

2.5. GO and KEGG enrichment analysis of DEGs

The analysis of Gene Ontology (GO) enhances our understanding of the biological functions associated with differentially expressed genes by identifying significant enrichments compared to the genomic background. Initially, we mapped all differentially expressed genes to terms in the GO database and calculated their respective gene counts [26]. We then observed a significant enrichment of genes when compared to the entire genome background. To identify significantly enriched terms (Padj≤0.05), clusterProfiler’s enricher function employs a hypergeometric distribution test as part of its GO enrichment analysis method.

Pathway enrichment analysis allows for the identification of crucial biochemical metabolic pathways and signal transduction pathways associated with differentially expressed genes. KEGG (Kyoto Encyclopedia of Genes and Genomes) serves as the primary public database for pathway information, integrating data on genome chemistry and system function. The unit of analysis in pathway significant enrichment analysis is the KEGG Pathway [27]. A hypergeometric test was used to identify pathways where differentially expressed genes were significantly enriched compared to all annotated genes. Similarly, a padj threshold less than 0.05 was employed to determine significant enrichment in KEGG pathway analysis.

2.6. Analysis of WGCNA

The Weighted Gene Co-expression Network Analysis (WGCNA) algorithm was used to define the correlation matrix of gene co-expression and the adjacency function formed by the gene network. Subsequently, various coefficients were calculated for different nodes, and a hierarchical clustering tree was constructed with distinct branches representing different gene modules [28]. Then, genes were identified based on their connectivity degree within each module, followed by an investigation into the relationship between these core genes and the effects of nano-selenium and organic selenium on leaves.

2.7. Quantitative real-time PCR

To ensure the reliability of RNA-Seq results, we validated nine differentially expressed genes (DEGs) identified by RNA-Seq through quantitative real-time polymerase chain reaction (qRT-PCR). We confirmed the expression patterns of these DEGs in leaves treated with organic selenium and nano-selenium by using SiActin as an internal reference and designing primers based on the sequences of DEGs. The RNA from the constructed library was reverse-transcribed into cDNA using PrimeScript TM 1st strand cDNA Synthesis Kit, followed by incorporating fluorescent dyes into the PCR reaction system. Real-time detection of fluorescence signals for each cycle product in PCR amplification generated a fluorescence amplification curve that underwent quantitative analysis [27]. Finally, gene expression levels were calculated and analyzed using the 2^-ΔΔCt method [29].

2.8. Statistical analysis

The data were sorted, and graphs were generated using Excel 2019 software. Statistical analysis was conducted utilizing SPSS 16.0 software (SPSS Inc., Chicago, IL, USA). The significant differences between the groups were determined using one-way ANOVA, followed by Duncan’s multiple range test.

3. Results

3.1. Photosynthetic pigments and photosynthetic indexes of Suzhouqing leaves

The results presented in Table 2 demonstrate that DN5, DN10, and DN20 increased chlorophyll a content by 68.97%, 41.38%, and 21.84% respectively compared to CK. Additionally, DO5, DO10, and DO20 resulted in an increase of chlorophyll a content by 1.15%, 20.69%, and 6.90% respectively. Furthermore, DN5, DN10, DN20, DO5, and DO10 led to an increase in chlorophyll b content by 46.65%, 23.91%, 6.52%, 2.17%, and 13.04% respectively; however, there was a decrease of 8.70% observed with the application of DO20 treatment. Moreover, DN5 and DN10 showed an increase in carotenoid content observed with DN20, DO5, and DO10 treatments. Conversely, a significant decrease of 39.38% was noted with the application of DO20.

As shown in Table 3, the stomatal conductance (Gs) of leaves initially increased and subsequently decreased with the increase in nano-selenium concentration, yet they remained consistently higher than those of CK, exhibiting an increase of 59.51%, 4.03%, and 4.08% respectively. The impact of nano-selenium on intercellular CO₂ concentration (Ci) was consistent with its effect on Gs. These findings suggest that selenium application enhances stomatal conductance, reduces stomatal resistance, and elevates intercellular carbon dioxide concentration in leaves. In comparison to CK, the net photosynthetic rate (Pn) of leaves exhibited a significant increase by 50.37%, 106.83%, and 108.11% respectively, while the transpiration rate (Tr) experienced a rise of 37.44% (DN5), followed by decreases of 2.10% (DN10), and 2.43% (DN20). Notably, these differences reached their maximum value during DN5 treatment.

3.2. Antioxidant properties of leaves

Figure 1 shows that he SOD activity was reduced by spraying nano-selenium, with DN5, DN10, and DN20 decreasing by 20.91%, 62.57%, and 29.68% respectively. Similarly, the SOD activity decreased when organic selenium was sprayed, with DO5, DO10, and DO20 decreasing by 17.47%, 53.36%, and 43.82% respectively. The POD activity increased by 14.42% and 35.79% under DN5 and DN20 respectively while it decreased by 31.42% under DN10. Compared to CK, the three groups treated with organic selenium spray showed a wider range of changes: DO5 and DO20 decreased by 23.99% respectively while DO10 increased by 10.43%. The CAT activity of the three groups sprayed with nano-selenium was reduced by 42.66%, 66.44% and 71.22% respectively compared to CK. Additionally, there was a decrease in CAT activity of 8.92%, 69.36% and 64.60% with organic selenium. The MDA content showed significant increases of 22.33%, 5.63%, and 7.71% with DN5, DN10, and DN20 treatments, respectively. Additionally, the application of organic selenium resulted in a wide range of changes: it decreased by 20.78% under DO5 treatment but increased by 43.03% and 0.81% under DO10 and DO20 treatments, respectively. Spraying nano-selenium increased the GSH content by percentages of 23.82%, 28.00%, and 4.98%, respectively compared to CK. The content of GSH also increased with the application of organic selenium by percentages of 34.89%, 12.77% and 9.55%,respectively compared to CK. Spraying low concentration and medium concentration (DN5 and DN10) nano-selenium decreased GSH content by 37.87% and 10.49%, respectively while high concentration (DN20) increased it by 63.77%. However, spraying organic selenium resulted in an increase in GSH content by 26.78%, 3.91%, and 0.73%, respectively.

3.3. Weight and total selenium content of Suzhouqing leaves

The results presented in Table 4 demonstrate that all treatments significantly increased the weight and total selenium content of Suzhouqing leaves. Notably, the DN10 treatment exhibited the highest values, with 67.30 g and a total selenium content of 42.37 μg kg⁻¹. These values were respectively 1.38 times and 4.15 times higher than those observed in the control group.

3.4. Nutritional quality of Suzhouqing leaves

As shown in Table 5, the soluble sugar content was increased by 70.43%, 9.43%, and 76.72% in comparison with CK, while it decreased by 49.94%, 10.28%, and increased 43.16% when treated with organic selenium. The soluble protein content in DN5 and DN10 increased by ratios of 54.61% and 3.51%. Spraying organic selenium resulted in increase of soluble protein content by 9.10% and 38.88% (DO5 and DO20). A low concentration of nano-selenium (DN5) led to increase in amino acid content by 7.51%; however, DN10 and DN20 increased it by 20.81% and 98.84%. DO5 and DO10 caused decreases of 6.36% and increased 5.78% respectively, whereas DO20 increased it by 16.76%. Treatment with DN5 increased the flavonoid content by 100.00%; however, both DN10 and DN20 increased it by 27.42% and 1.61%. After spraying organic selenium, the flavonoid content decreased by 8.06%, 11.29% and 70.97%, respectively.

3.5. Analysis of the transcriptomic changes of nano-selenium and organic selenium in response to drought stress

3.5.1. Sequencing results and quality analysis

Transcriptome sequencing has been demonstrated in numerous species as a highly effective strategy for unraveling their global expression networks under stress conditions. To investigate the molecular mechanisms associated with nano-selenium and organic selenium in response to drought stress, transcriptome sequencing was conducted on collected samples, and transcriptome libraries were constructed using total RNA extracted from leaves of both the control group and treatment groups. The transcriptome libraries underwent Illumina sequencing, with treatments including a control group, organic selenium and nano-selenium. Subsequently, data filtering was performed on the raw data resulting in 297,176,114 clean sequences obtained for analysis. Each sample had a sequencing volume exceeding 5.50 G. Notably, over 99.80% of the total base count consisted of Q30 bases as presented in Table 6.

3.5.2. Correlation examination between samples

The correlation of gene expression levels between samples serves as a crucial indicator for assessing the experiment’s reliability and the rationality of sample selection. A higher correlation coefficient, closer to 1, signifies greater similarity in expression patterns among the samples. Intra-group and inter-group correlation coefficients were computed based on TPM values of all genes within each sample, followed by generating a heat map to visually represent inter-group differences and intra-group duplications. With a minimum squared Pearson correlation coefficient value of 0.85, it is evident that gene expression levels across replicates within each group exhibit high consistency.

3.5.3. Differential expressed genes (DEGs) statistics

The differential significance analysis was conducted for each comparison combination, resulting in the identification of a list of differentially expressed genes (DEGs) based on the criteria of |log₂(fold change)| ≥ 1 and adjusted p-value ≤ 0.05. This list is presented in the table below. A total of 58,852 DEGs were identified between DN5 and DCK, with 21,358 up-regulated DEGs (log₂(fold change) ≥ 1) and 21,530 down-regulated DEGs (log₂(fold change) ≤ -1).

Table 7. List of DEGs.

Group	All	Up-regulated	Down-regulated
DN5vsDCK	58,852	21,358	21,530
DN10vsDCK	60,711	22,397	23,018
DN20vsDCK	58,566	21,519	22,818
DO5vsDCK	57,070	20,651	24,828
DO10vsDCK	57,473	20,059	22,953
DO20vsDCK	59,529	21,324	21,444

Table 7. List of DEGs.

Group	All	Up-regulated	Down-regulated
DN5vsDCK	58,852	21,358	21,530
DN10vsDCK	60,711	22,397	23,018
DN20vsDCK	58,566	21,519	22,818
DO5vsDCK	57,070	20,651	24,828
DO10vsDCK	57,473	20,059	22,953
DO20vsDCK	59,529	21,324	21,444

3.5.4. Volcano map of DEGs

The differential expression analysis focuses on identifying genes that are differentially expressed (DEGs) between samples and subsequently functionally analyzing these genes. Based on the criteria of |log₂(fold change)|≥1 and FDR<0.05 for significant DEGs, our results revealed a total of 58,852 and 60,711 and 58,566 significantly differentially expressed genes between DN5, DN10, and DN20 compared to DCK respectively (Figure 2a–c). Additionally, we observed significant differential expression in 57,070 and 57,473 and 59,529 genes between DO5, DO10 and DO20 compared to DCK respectively (Figure 2d–f).

3.5.5. Venn diagram of DEGs

The Venn diagram illustrates the intersection of differentially expressed genes across various comparison combinations, enabling the identification of common or unique differential genes within specific comparisons. The 36,078 DEGs that are common exhibit consistent expression across all combinations, as illustrated in Figure 3.

3.5.6. GO annotation and enrichment analysis of differentially expressed genes

The functional annotation of DEGs between the nano-selenium and organic selenium treatment groups was conducted. Subsequently, GO functional enrichment analysis was separately performed for the up-regulated and down-regulated DEGs. Statistical tests were carried out using corrected P-values, and a histogram depicting the distribution of GO terms was generated (Figure 4).

In terms of biological processes, it was observed that up-regulated genes tended to be involved in photosynthesis-related processes such as light harvesting (GO: 0009765), photosynthetic electron transport chain (GO: 0009767), green metabolic system process (GO: 0015994), cellular divalent inorganic cation homeostasis (GO:0072503), green metabolic process (GO: 0072503), and regulation of chlorophyll biosynthetic process (GO:0010380). On the other hand, down-regulated genes were primarily associated with regulation of protein localization to nucleus (GO: 1900180) and regulation of dephosphorylation events (GO: 0035303). Regarding cellular components, major up-regulated genes were found to be related to chloroplast thylakoid membrane protein complex formation (GO: 0098807), chloroplast thylakoid lumen organization (GO: 0009543), transcriptionally active chromatin maintenance (GO: 0035327); whereas down-regulated genes mainly belonged to cajal body organization (GO: 0015030) and tubulin complex assembly (GO: 0045298). In terms of molecular functions, significant expression changes were observed for pigment binding activity (GO: 0031409), structural constituent of cytoskeleton (GO: 000520) and poly(U) RNA binding activity (GO: 000826). On the contrary, down-regulated gene sets mainly included snoRNA binding activity (GO: 003051) and acid-thiol ligase activity (GO: 0016878).

Figure 4. GO enrichment histogram of differential genes. Note: B: Biological process; C: Cellular component; M: Molecular function; B1: Cell plate assembly; B2: Cell plate formation involved in plant-type cell wall biogenesis; B3: Cellular divalent inorganic cation homeostasis; B4: Cellular response to gravity; B5: Cellular zinc ion homeostasis; B6: Central nervous system neuron development; B7: Chlorophyll metabolic process; B8: Circadian regulation of calcium ion oscillation; B9: Circadian regulation of gene expression; B10: Clathrin-dependent endocytosis; B11: Defense response to Gram-negative bacterium; B12: Ethanolamine-containing compound metabolic process; B13: Long-day photoperiodism; B14: Negative regulation of GTPase activity; B15: Nitric oxide biosynthetic process; B16: Nitric oxide metabolic process; B17: NLS-bearing protein import into nucleus; B18: Nucleotide-excision repair, DNA gap filling; B19: Photosynthesis, light harvesting; B20: Photosynthesis, light harvesting in photosystem I; B21: Photosynthetic electron transport chain; B22: Photosynthetic electron transport in photosystem II; B23: Photosystem I stabilization; B24: Photosystem II repair; B25: Plastid translation; B26: Positive regulation of protein import; B27: Positive regulation of protein import into nucleus; B28: Positive regulation of protein localization to nucleus; B29: Primary miRNA processing; B30: Protein repair; B31: Red, far-red light phototransduction; B32: Regulation of chlorophyll biosynthetic process; B33: Regulation of dephosphorylation; B34: Regulation of DNA-templated transcription in response to stress; B35: Regulation of photosynthesis, light reaction; B36: Regulation of protein dephosphorylation; B37: Regulation of protein import; B38: Regulation of protein import into nucleus; B39: Regulation of protein localization to nucleus; B40: Regulation of tetrapyrrole biosynthetic process; B41: Response to low fluence blue light stimulus by blue low-fluence system; B42: Response to low light intensity stimulus; B43: Response to strigolactone; B44: Ribosomal protein import into nucleus; B45: Salicylic acid catabolic process; B46: Spindle assembly involved in meiosis; C1: Cajal body; C2: Chloroplast thylakoid lumen; C3: Chloroplast thylakoid membrane protein complex; C4: Light-harvesting complex; C5: Nuclear dicing body; C6: Photosystem; C7: Plastid ribosome; C8: Plastid thylakoid lumen; C9: Protein phosphatase type 2A complex; C10: Proton-transporting V-type ATPase, V0 domain; C11: thylakoid light-harvesting complex; C12: Transcriptionally active chromatin; C13: Tubulin complex; C14: Vacuolar proton-transporting V-type ATPase, V0 domain; M1: Acid-thiol ligase activity; M2: Alcohol binding; M3: Chlorophyll binding; M4: Glutaminase activity; M5: Isoprenoid binding; M6: MAP-kinase scaffold activity; M7: Methionine adenosyltransferase activity; M8: Pigment binding; M9: Poly(U) RNA binding; M10: Pyrophosphate hydrolysis-driven proton transmembrane transporter activity; M11: Signaling adaptor activity; M12: SnoRNA binding; M13: Structural constituent of cytoskeleton.

Figure 4. GO enrichment histogram of differential genes. Note: B: Biological process; C: Cellular component; M: Molecular function; B1: Cell plate assembly; B2: Cell plate formation involved in plant-type cell wall biogenesis; B3: Cellular divalent inorganic cation homeostasis; B4: Cellular response to gravity; B5: Cellular zinc ion homeostasis; B6: Central nervous system neuron development; B7: Chlorophyll metabolic process; B8: Circadian regulation of calcium ion oscillation; B9: Circadian regulation of gene expression; B10: Clathrin-dependent endocytosis; B11: Defense response to Gram-negative bacterium; B12: Ethanolamine-containing compound metabolic process; B13: Long-day photoperiodism; B14: Negative regulation of GTPase activity; B15: Nitric oxide biosynthetic process; B16: Nitric oxide metabolic process; B17: NLS-bearing protein import into nucleus; B18: Nucleotide-excision repair, DNA gap filling; B19: Photosynthesis, light harvesting; B20: Photosynthesis, light harvesting in photosystem I; B21: Photosynthetic electron transport chain; B22: Photosynthetic electron transport in photosystem II; B23: Photosystem I stabilization; B24: Photosystem II repair; B25: Plastid translation; B26: Positive regulation of protein import; B27: Positive regulation of protein import into nucleus; B28: Positive regulation of protein localization to nucleus; B29: Primary miRNA processing; B30: Protein repair; B31: Red, far-red light phototransduction; B32: Regulation of chlorophyll biosynthetic process; B33: Regulation of dephosphorylation; B34: Regulation of DNA-templated transcription in response to stress; B35: Regulation of photosynthesis, light reaction; B36: Regulation of protein dephosphorylation; B37: Regulation of protein import; B38: Regulation of protein import into nucleus; B39: Regulation of protein localization to nucleus; B40: Regulation of tetrapyrrole biosynthetic process; B41: Response to low fluence blue light stimulus by blue low-fluence system; B42: Response to low light intensity stimulus; B43: Response to strigolactone; B44: Ribosomal protein import into nucleus; B45: Salicylic acid catabolic process; B46: Spindle assembly involved in meiosis; C1: Cajal body; C2: Chloroplast thylakoid lumen; C3: Chloroplast thylakoid membrane protein complex; C4: Light-harvesting complex; C5: Nuclear dicing body; C6: Photosystem; C7: Plastid ribosome; C8: Plastid thylakoid lumen; C9: Protein phosphatase type 2A complex; C10: Proton-transporting V-type ATPase, V0 domain; C11: thylakoid light-harvesting complex; C12: Transcriptionally active chromatin; C13: Tubulin complex; C14: Vacuolar proton-transporting V-type ATPase, V0 domain; M1: Acid-thiol ligase activity; M2: Alcohol binding; M3: Chlorophyll binding; M4: Glutaminase activity; M5: Isoprenoid binding; M6: MAP-kinase scaffold activity; M7: Methionine adenosyltransferase activity; M8: Pigment binding; M9: Poly(U) RNA binding; M10: Pyrophosphate hydrolysis-driven proton transmembrane transporter activity; M11: Signaling adaptor activity; M12: SnoRNA binding; M13: Structural constituent of cytoskeleton.

3.5.7. KEGG annotation and enrichment analysis of differentially expressed genes

(1) KEGG pathway: photosynthesis-antenna proteins

Pathway-based analysis is crucial for gaining a deeper understanding of the biological functions of genes. KEGG serves as a prominent public database for pathway information, and hypergeometric analysis is employed to identify significant enrichment of pathways.

Usually, PS I light trapping pigment protein (LHC) contains Lhca1, Lhca2, Lhca3 and Lhca4 proteins. PS II light-trapping pigment protein (LHCII) contains Lhcb1, Lhcb2, Lhcb3, Lhcb4, Lhcb5 and other proteins, and the molecular weight of light-trapping pigment protein is between 20-29 kD. There are two kinds of chlorophyll a/b pigment protein complexes in the inner capsule membrane, the chlorophyll a/b pigment protein complex (LHCI) of the photocatcher antenna and the chlorophyll a/b pigment protein complex of the photocatcher antenna (LHCII) of the photocatcher antenna (PS II). Since LHCII combines 50% of the total leaf chlorophyll, it plays an important role in the absorption of light energy in photosynthesis. The process of plant photosynthesis is highly susceptible to the effects of drought stress, which can impede photochemical activity and disrupt the equilibrium between light energy absorption and utilization. PSII and PSI serve as the reaction center complexes that drive the photocooperative photoreaction. The biogenesis of PSII involves the co-assembly of at least 20 distinct peptides and a multitude of cofactors, while PSI is a multi-subunit protein complex located in the chloroplast thylakoid membrane. Under drought stress conditions, damage to photoinhibition relies on the rate of de novo assembly of all subunits and cofactors in PSII, indicating that downregulation of photosynthetic genes may lead to decreased photosynthesis.

The light-harvesting complex (LHC) serves as a photoreceptor, efficiently capturing and transferring excitation energy to closely associated photosystems. When the q-value is equal to or less than 0.05, there are notable differences observed in the drought response of photosynthetic antenna (light pigment protein complexes) and their impact on light-harvesting protein complex (LHC) and photosynthetic capacity. Figure 5, which refers to schematic diagram 4 from Kanehisa laboratories depicting photosynthetic antenna proteins, highlights the significantly differentially expressed genes, while Table 8 provides a list of these corresponding differentially expressed genes. There was a notable increase in the abundance of phototrapping chlorophyll a/b pigment antenna proteins in PS II (Table 8), indicating that both nano-selenium and organic selenium treatments augmented leaf photoelectron trapping capacity, thus significantly enhanced the efficiency of photosystem I (PS I) and photosystem II (PS II) compared to the control group.

(2) KEGG pathway: thiamine metabolism

Thiamine is a bicyclic compound consisting of thiazole and pyrimidine moieties. In plants, thiamine exists primarily in three forms: free thiamine, thiamine monophosphate and thiamine pyrophosphate. During the biosynthesis process of thiamine, the pyrimidine and thiazole components are synthesized independently coupled to form thiamine pyrophosphate. Pyrimidine synthetase (THIC), thiazole synthetase (THI1), thiamine phosphate synthetase (TH1), and thiamine pyrophosphokinase (TPK) serve as key enzymes in the pathway regulating plant’s synthesis of this essential vitamin (Figure 6). Thiamine plays a crucial role in plant response to drought stress, and thiamine metabolism is shown in the figure below. As shown in Table 9, the expression of cysteine desulfurase, cysteine-dependent adenosine diphosphate thiazole synthase, adenylate kinase, and thiamine pyrophosphokinase was up-regulated upon treatment with organic selenium and nano-selenium.

(3) KEGG pathway: glutathione metabolism

The glutathione metabolic pathway plays a significant role in adapting to drought conditions, particularly in maintaining water homeostasis and reducing water sensitivity in leaves. This has been elucidated through key pathways such as reactive oxygen species (ROS) clearance and signal transduction, as well as the up-regulation of antioxidant enzyme gene expression (Figure 7) and Table 10.

(4) KEGG pathway: selenocompound metabolism

We used bioinformatics analysis to explore the leaf genome and determine the comprehensive metabolic pathway of selenium compounds (Figure 8). The metabolism of selenocompounds begins with the conversion of selenocysteine through cystathionine gamma-synthase (CGS), resulting in the release of seleno-cystathionine. The subsequent compound undergoes metabolic transformation via cysteine-S-conjugate beta-lyase (CBL), leading to the liberation of selenohomocysteine. This resulting compound reacts with 5-methyltetrahydrofolate-homocysteine methyltransferase (MET), resulting in the liberation of seleno-methionine. Selenomethionine can be further metabolized into either seleno-methionyl-tRNA or methyl-selenol. Alternatively, methyl-selenol can also be formed from the reaction between methyl-selenic acid and thioredoxin reductase or Se-methyl-seleniumcysteine through cystathionin gamma-lyase (CTH). Compared to the control group, there were significant differences in gene expression between nano-selenium treatment and organic selenium treatment, as shown in Table 11.

3.7. WGCNA of DEGs and identification of hub gene

The WGCNA assumes a scale-free gene network, defines the correlation matrix and adjacency function for gene co-expression, and calculates different coefficients for various nodes. Subsequently, a hierarchical clustering tree is constructed. Different branches of the cluster tree represent distinct gene modules characterized by high intra-module gene co-expression and varying scores indicating low inter-module gene co-expression. The correlation between the module and the sample is shown in Figure 9. The results revealed that DN5 exhibited the highest correlation (0.88) with the MEsalmon module, while DN10 displayed the highest correlation (0.93) with the MEtan module. Additionally, DN20 demonstrated the strongest association (0.96) with the MEdarkred module, and DO5 showed a significant (0.95) with the MEmagenta module. Moreover, DO10 exhibited the highest correlation (0.95) with the MEbrown module, whereas DO20 displayed a strong association (0.96) with the MEdarkorange module.

The expression trend of genes in the module across different samples is illustrated in Figure 10, and the expression patterns of the aforementioned six module eigenvalues are analyzed. Subsequently, the connectivity degree (commonly referred to as connectivity or degree, denoted by k) between a gene within the module and other genes is calculated. Generally, genes with high connectivity rankings (k values) within a module can be considered as hub genes due to their central position. The results of connectivity are presented in the Table 12, and a total of nine hub genes were obtained: LOC103837231, Adenylate kinase 2; LOC103858064, glutathione peroxidase 3; LOC103842805, glutathione S-transferase U25; LOC103830524, glutathione S-transferase U20; LOC103844552, thiamine pyrophosphokinase 1; LOC103872309, ABC transporter G family member 35; LOC103867761, ABC transporter G family member 1; LOC103844798, ATP-dependent zinc metalloprotease FTSH 11; LOC103866104, cysteine synthase.

The observed gene expression changes by WGCNA were consistent with the results of up-regulated genes in the KEGG pathway of DEGs as above by transcriptome sequencing analysis, thereby validating the accuracy and reliability of the transcriptome sequencing data. Additionally, differential expression of these genes was observed between the experimental and control groups. Consequently, these genes can be considered crucial selenium response-related genes that play a pivotal role in the molecular response mechanism of leaf towards exogenous selenium treatment.

Table 12. Gene connectivity results for the six modules.

Gene id	Gene name	Function	kTotal	kWithin	kOut	kDiff	Module
Cluster-1111.26333	LOC103837231	Adenylate kinase 2	128.913278	53.85615	75.05713	21.201	Salmon
Cluster-1111.33663	LOC103858064	Glutathione peroxidase 3	266.4526	169.3169	97.1357	72.18123	Salmon
Cluster-1111.62702	LOC103842805	Glutathione S-transferase U25	788.70	573.26	215.44	357.81	Magenta
Cluster-1111.32944	LOC103830524	Gglutathione S-transferase U20	216.4186	126.5464	89.87223	36.67415	Salmon
Cluster-1111.94080	LOC103844552	Thiamine pyrophosphokinase 1	77.49846	14.04587	63.45259	49.4067	Darkorange
Cluster-1111.21138	LOC103872309	ABC transporter G family member 35	470.7212	281.3282	189.393	91.93526	Darkred
Cluster-1111.51912	LOC103867761	ABC transporter G family member 1	317.63	200.26	117.38	82.88	Darkorange
Cluster-1111.47151	LOC103844798	ATP-dependent zinc metalloprotease FTSH 11	741.37	564.38	176.99	387.39	Tan
Cluster-1111.43447	LOC103866104	Cysteine synthase	83.49206	27.36962	56.12244	28.7528	Brown4

Note: kTotal: total gene connectivity; kWithin: the connectivity of genes within modules; kOut: the connectivity of genes outside the module, kTotal minus kWithin; kDiff: The difference between kWithin and kOut.

Table 12. Gene connectivity results for the six modules.

Gene id	Gene name	Function	kTotal	kWithin	kOut	kDiff	Module
Cluster-1111.26333	LOC103837231	Adenylate kinase 2	128.913278	53.85615	75.05713	21.201	Salmon
Cluster-1111.33663	LOC103858064	Glutathione peroxidase 3	266.4526	169.3169	97.1357	72.18123	Salmon
Cluster-1111.62702	LOC103842805	Glutathione S-transferase U25	788.70	573.26	215.44	357.81	Magenta
Cluster-1111.32944	LOC103830524	Gglutathione S-transferase U20	216.4186	126.5464	89.87223	36.67415	Salmon
Cluster-1111.94080	LOC103844552	Thiamine pyrophosphokinase 1	77.49846	14.04587	63.45259	49.4067	Darkorange
Cluster-1111.21138	LOC103872309	ABC transporter G family member 35	470.7212	281.3282	189.393	91.93526	Darkred
Cluster-1111.51912	LOC103867761	ABC transporter G family member 1	317.63	200.26	117.38	82.88	Darkorange
Cluster-1111.47151	LOC103844798	ATP-dependent zinc metalloprotease FTSH 11	741.37	564.38	176.99	387.39	Tan
Cluster-1111.43447	LOC103866104	Cysteine synthase	83.49206	27.36962	56.12244	28.7528	Brown4

Note: kTotal: total gene connectivity; kWithin: the connectivity of genes within modules; kOut: the connectivity of genes outside the module, kTotal minus kWithin; kDiff: The difference between kWithin and kOut.

3.8. Fluorescence quantitative PCR analysis

Nine genes were selected from the differentially expressed genes in the transcriptome for validation analysis using real-time fluorescence quantitative PCR, and the correlation coefficient was 0.95 (calculated by 2^-ΔΔCt value). The results of qRT-PCR verification exhibited a significant positive correlation with RNA-Seq results at a level of p< 0.01, indicating that the test outcomes were reliable (Figure 10).

Figure 10. qRT-PCR identification and RNA-Seq analysis. Note: A: LOC103837231, adenylate kinase 2; B: LOC103858064, glutathione peroxidase 3; C: LOC103842805, glutathione S-transferase U25; D: LOC103830524, glutathione S-transferase U20; E: LOC103844552, thiamine pyrophosphokinase 1; F: LOC103872309, ABC transporter G family member 35; G: LOC103867761, ABC transporter G family member 1; H: LOC103844798, ATP-dependent zinc metalloprotease FTSH 11; I: LOC103866104, cysteine synthase.

Figure 10. qRT-PCR identification and RNA-Seq analysis. Note: A: LOC103837231, adenylate kinase 2; B: LOC103858064, glutathione peroxidase 3; C: LOC103842805, glutathione S-transferase U25; D: LOC103830524, glutathione S-transferase U20; E: LOC103844552, thiamine pyrophosphokinase 1; F: LOC103872309, ABC transporter G family member 35; G: LOC103867761, ABC transporter G family member 1; H: LOC103844798, ATP-dependent zinc metalloprotease FTSH 11; I: LOC103866104, cysteine synthase.

4. Discussion

4.1. Effects of organic selenium and nano-selenium on leaf photosynthesis under drought stress

Chlorophyll is the main pigment that maintains plant photosynthesis, directly affecting the strength of plant photosynthetic capacity. Yao et al. [30] found that exogenous selenium could increase chlorophyll content in wheat seedling leaves under drought stress. Chen et al. [31] conducted experiments on tobacco by controlling water and simulating drought stress, finding that exogenous selenium could promote the synthesis of photosynthetic pigments in tobacco leaves under drought stress and increase chlorophyll content. Yue [32] discovered that exogenous selenium could alleviate chlorophyll destruction in N51 seedling leaves under drought stress, better helping corn seedlings resist damage caused by drought stress and improving their drought resistance. The results of this study are essentially consistent with those of the aforementioned scholars. However, in this paper, digital gene expression profile technology was employed to analyze differential gene expression. The findings revealed a significant increase in the expression levels of photosynthetic antenna protein regulatory genes and major regulatory genes of photosynthesis upon spraying nano-selenium and organoselenium, leading to an improvement in photosynthesis. Additionally, exogenous selenium treatment resulted in a notable increase in the net photosynthetic rate, transpiration rate, and stomatal conductance of leaves, indicating an enhancement in photosynthetic capacity that aligns with the results obtained from gene expression profile analysis.

4.2. Effects of organic selenium and nano-selenium on regulation of antioxidant system of leaves under drought stress

The application of selenium has been shown to increase the activity of antioxidant enzymes and reduce oxidative damage in wheat under water stress [33]. Fan [34] demonstrated that selenium supplementation significantly decreased the accumulation of H₂O₂, O₂^-, and MDA in plant leaves and roots under drought and salt stress, which is consistent with the findings of Rady et al. [35]. Yue [32] found exogenous selenium could significantly enhance the activity of antioxidant enzymes in N51 seedling leaves and roots under drought stress, reduce reactive oxygen species levels in vivo, maintain normal tissue physiological activities, and improve seedling drought resistance.

Ahmad et al. [36] discovered that foliar application of selenium and priming with selenium significantly enhanced the activities of SOD, POD, and CAT in oilseed crops Camelina (Camelina sativa L.) and Canola (Brassica napus L.) under drought stress conditions. Sardari et al. [37] reported that selenium and nano-selenium acted as stimulators for antioxidant enzyme activities in wheat plants, thereby enhancing plant tolerance to drought stress. Zahedi et al. [38] demonstrated that foliar spraying of 10 nm and 50 nm Se-NPs alleviated the detrimental effects of drought on pomegranate leaves and fruits by reducing stress-induced lipid peroxidation and H₂O₂ content, while enhancing the activity of antioxidant enzymes. Moloi and Khoza [39] observed that selenium induced an increase in ascorbate peroxidase activity during the flowering stage in edamame plants under drought stress.

The findings of this study are consistent with previous research: the application of organic selenium and nano-selenium resulted in an increase in SOD activity, POD activity, and GSH content. Additionally, transcriptome sequencing revealed that this could be attributed to selenium acting as a cofactor for glutathione peroxidase (GPX) and participating in the antioxidant defense system. The up-regulation of glutathione S-transferase and glutathione oxidase genes reduced oxidative damage in plants under stress, indirectly reducing the need for antioxidant enzymes to eliminate ROS.

Based on research findings on thiamine metabolism, we hypothesize that the application of nano-selenium and organic selenium can help establish a balance between the production and elimination of reactive oxygen species by thiamine. The clearance of reactive oxygen species may be associated with direct coupling of its molecules, while its metabolite thiamine pyrophosphate is stored as a coenzyme to enhance mitochondrial oxidation state under foreign stress stimuli. This enables cells to release reactive oxygen species more efficiently when stimulated by external factors, thereby inducing downstream signaling molecules for plant resistance. However, further investigation is needed to explore the intricate relationship between thiamine and mitochondrial oxidative state.

4.3. Effects of organic selenium and nano-selenium on weight, total selenium content and leaf nutritional quality under drought stress

The metabolism and regulatory mechanism of selenium in plants constitute a complex gene expression network. Various intermediates synthesized in the pathway of selenium metabolism may participate in regulating plant growth and development. In this study, we found that nano-selenium and organic selenium could enhance the accumulation of soluble sugars, soluble proteins, free amino acids, and ascorbic acid in leaves. Soluble sugar acts as a cell osmotic regulator while ascorbic acid serves as a potent antioxidant, both playing crucial roles in maintaining plant tolerance to environmental stress. Soluble proteins and free amino acids are vital constituents of numerous enzymes involved in diverse metabolic processes. These substances represent important primary metabolites within plants and are frequently employed as representative indicators for evaluating the growth status and nutritional quality of horticultural products. These compounds can serve as osmoregulatory factors to maintain osmotic pressure balance in plants and protect plant cells in stressful environments [40].

In recent years, studies have demonstrated that an adequate supply of selenium can enhance crop yield and increase grain selenium content. This finding is consistent with the results of our study, which revealed that both nano-selenium and organic selenium treatments promote the growth and development of Suzhouqing. These findings align with previous research conducted by Li et al. [41], Guo et al. [42], Zhang et al. [43], Gao et al. [44], Kong [29] on millet, as well as Song et al. [45] and Peng et al. [46] on wheat cultivation practices.

Transcriptome analysis of drought-stressed cotton revealed that differentially expressed genes were significantly enriched in the osmoregulatory system [47]. Transcriptome sequencing was employed to elucidate the absorption, transport, and metabolism of selenium in Cymbidium japonicum. The investigation identified the relevant pathways and key enzyme genes involved in selenium metabolism while exploring the mechanism underlying selenium enrichment. Hu et al. [48] discovered 16 crucial genes belonging to the ABC transporter family in tea plants, including 4 members from the ABCB subfamily, 4 from the ABCC subfamily, 5 from the ABCG subfamily, as well as 2 members from both ABCA and ABCF subfamilies. It is hypothesized that ABC transporters may participate in selenide transportation within tea plants. Kong [29] examined the physiological and biochemical effects of applying a biological milkshake on yield and found that four out of five differential genes related to ABC transporters were up-regulated while one was down-regulated. In our study, multiple key genes associated with ABC transporters were primarily identified within the ABCG subfamily; these differential genes exhibited up-regulation without any down-regulation. This finding contrasts with Kong’s results [29], suggesting potential tissue specificity differences regarding ABC transporter function across different species.

In this study, we observed varying degrees of up-regulation in the key genes controlling the expression of cystathionine gamma-lyase and 5-methyltetrahydrofolate-homocysteine methyltransferase. These findings suggest that the leaf enhances selenium absorption by up-regulating the expression of relevant genes in response to selenium supply in its environment. Furthermore, our results are consistent with previous studies conducted by Freeman [49] and Hoewyk [50], which also identified consistent patterns in the hub genes involved in cysteine synthase. However, further research is required to investigate how these genes directly participates in the absorption, transport and metabolism of selenium.

5. Conclusions

The underlying physiological and molecular mechanisms of organic selenium and nano-selenium on Suzhouqing were examined using a transcriptomics approach, and the following conclusions were drawn. Under drought conditions, the beneficial effects of organic selenium and nano-selenium on the weight and quality of ‘Suzhouqing’ were primarily manifested in the following four aspects. Firstly, there was an enhancement in photosynthetic capacity through upregulation of light-trapping pigment proteins such as Lhca1, Lhca2, Lhca3, Lhca4, Lhcb1, Lhcb2, Lhcb3, Lhcb4 and Lhcb5. Secondly, effective regulation of reactive oxygen species (ROS) homeostasis was achieved by activating the antioxidant system via up-regulation of glutathione S-transferase. Thirdly, water homeostasis was maintained through glutathione oxidase activity. Lastly, increased expression levels of ABC transporter, adenylate kinase and cysteine desulphurase contributed to elevated total selenium content.

The determination of the application effect of selenium on Suzhouqing under drought stress is not solely influenced by unilateral factors. Therefore, we will use multi-omics analysis to further investigate how Suzhouqing responds to drought stress when selenium is applied, ensuring its adaptation to a drought environment while simultaneously enhancing weight.