Effects of Foliar Spraying of Organic Selenium and Nano-Selenium Fertilizer on Pak Choi ( Brassica chinensis var. pekinensis . cv. ‘Suzhouqing’) under Low Temperature Stress

: The effects of foliar spraying of organic selenium and nano-selenium fertilizer on pak choi ( Brassica chinensis var. pekinensis . cv. ‘Suzhouqing’) under low temperature were investigated. The impacts on plant growth, antioxidant capacities, and nutritional proﬁles were studied. Exogenous selenium was applied at three rates: 5, 10, 20 mg L − 1 , and RNA-Seq technology was used to sequence the transcriptome of leaves. The study revealed that selenium inﬂuenced leaf weight and total selenium content through three main mechanisms. First, it protected photosynthetic pigments and boosted photosynthetic capacity by up-regulating LHca2, LHcb1, LHca1, and LHcb4. Second, it enhanced antioxidant capacity by elevating the expression of genes such as superoxide dismutase and monodehydroascorbic acid. Third, it facilitated selenium absorption through endocytosis, transported selenium using the ABC transport gene family, and regulated selenium metabolism-related genes like cysteine synthetase and glutaredoxin. Nine hub genes identiﬁed with a weighted gene co-expression network analysis were closely associated with these mechanisms. The results of a functional enrichment analysis were consistent with those of a Gene Ontology enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway analysis conducted on DEGs, thus conﬁrming the reliability of these ﬁndings. Therefore, this study can provide scientiﬁc basis for pak choi production with selenium fortiﬁcation by selenium application.


Introduction
In recent years, low-temperature disasters caused by global climate change have occurred frequently, seriously affecting the normal growth of crops.Low temperature stress leads to a significant accumulation of reactive oxygen species in plants.And it subsequently results in membrane lipid peroxidation and the formation of a large number of harmful substances such as malondialdehyde [1].Pak choi (Brassica chinensis var.pekinensis cv.'Suzhouqing') is rich in fiber, vitamins, and mineral elements [2].It can be cultivated throughout the year, making it an essential vegetable supply [2].Although this variety exhibits stronger cold resistance compared to other Chinese cabbage varieties, its weight and quality are reduced under low-temperature conditions.
Selenium (Se), as a beneficial element, is a component of the antioxidant glutathione peroxidase (GSH-Px).Numerous studies have demonstrated that selenium can participate in the removal of lipid hydrogen, hydrogen peroxide, and other free radicals.The application of an appropriate selenium fertilizer can alleviate free radical-induced peroxidation damage and offset the negative effects of low temperatures on plant growth and development.Huang et al. investigated the physiological indices of strawberries under low temperature stress and found that selenium influenced the ascorbic acid (ASA)-glutathione (GSH) cycle to combat low temperature stress [3].Liu et al. discovered that selenium not only enhances cold resistance in tea plants but also improves tea quality by promoting the accumulation of tea polyphenols, amino acids, and sugars [4].Therefore, the objective of this study was to investigate the efficacy of applying organic selenium fertilizer and nano-selenium on Suzhouqing for enhancing its resistance against low temperature stress.
The application of exogenous selenium can not only induce plant growth and nutritional quality but also trigger the expression of a series of related genes in plants, thereby causing various physiological response mechanisms.With the advancements in molecular biology and bioinformatics, we have gained a deeper understanding of the gene regulatory networks involved in plant responses to selenium.In recent years, an increasing number of studies have utilized RNA sequencing technology to investigate the effects of selenium on plants.Rao et al. integrated sequencing results from PacBio SMRT and Illumina RNA-Seq platforms, analyzing a total of 948 differentially expressed genes, out of which 11 were found to be associated with selenium metabolism in the selenium-accumulating plant Cardamine violifolia [5].They revealed candidate genes and pathways involved in selenium metabolism.Zhou et al., after applying selenium treatment, sequenced the RNA of Cardamine hupingshanensis seedlings and identified 25 significantly responsive genes under selenium stress, showing significant enrichment [6].Guo et al. identified structural genes, phosphate transporters, and sulfate transporters that may be involved in selenium metabolism when analyzing the response of Pueraria lobata (Willd.)Ohwi to selenium stimulation at the molecular level [7].Cao et al. observed that differential genes in tea trees responding to selenium stimulation were significantly enriched in four pathways [8].However, the mechanism of Suzhouqing's response to nano-selenium and organic selenium under low temperature through RNA-Seq has not been reported so far.Therefore, we used transcriptome sequencing analysis to identify differentially expressed genes and analyze the metabolic pathways associated with these genes.The physiological, biochemical, and molecular regulatory mechanisms of leaf surface application of selenium under low temperature stress were finally revealed, providing a theoretical basis for further enhancing weight and quality while increasing selenium content.

Plant Materials and Experimental Conditions
The tested pak choi was Suzhouqing, purchased from Jiangxi Ruibao Seed Co., Ltd.(Jiujiang, China).Nano-selenium (a small red elemental nanoparticle of selenium with a particle diameter ranging from 20 to 60 nm) was developed by Zhongnong Selenium-rich Agricultural Technology Research Institute (Beijing, China), whereas organic selenium (mainly malt selenium, (CH 3 ) 2 Se) was developed by Northwest A&F University.
The experiment was conducted in the laboratory of Jiangxi Agricultural University.The seedlings were cultivated in November 2022.Once the seedlings reached the three-leaf stage, they were transplanted into planting pots (33 × 24.4 × 13.5 cm) on 18 November, with two plants per pot, using organotrophic soil as the substrate.The physical and chemical properties of the soil were as follows: organic matter content of 622.72 g kg −1 , alkali-hydrolyzed nitrogen content of 412.34 mg kg −1 , available phosphorus content of 102.20 mg kg −1 , available potassium content of 2.40 g kg −1 , total selenium content of 1.44 mg kg −1 , and a pH value of 5.7.Nitrogen, phosphorus, and potassium fertilizers were applied to all treatments at a rate of pure N at 3.50 g per pot and an N:P 2 O 5 :K 2 O ratio of 1:1.5:3 dissolved in water and applied once before seedling transplantation.After transplanting, watering twice a day during the seedling period was carried out to promote seedling growth.
The entire plant was sprayed with 20 mL of water, different concentrations of organic selenium, and nano-selenium solutions on 9 and 16 December 2022, respectively.The experimental setup was as follows: one control group (LCK: water) and six experimental groups consisting of nano-selenium at concentrations of 5 mg L −1 (LN5), 10 mg L −1 (LN10), and 20 mg L −1 (LN20), as well as organic selenium at concentrations of 5 mg L −1 (LO5), 10 mg L −1 (LO10), and 20 mg L −1 (LO20).Each treatment had three replicates.Subsequently, on December 19th, the seedlings were subjected to a low temperature stress treatment in an artificial climate chamber at a temperature of −2 • C, light intensity of 2000 lx, photoperiods alternating between 12 h/12 h, and relative humidity set at 75%.On 27 December, the plants were sprayed again with clean water and different concentrations of organic selenium and nano-selenium solutions.Finally, on 1 January 2023, the seedlings were removed from the artificial climate chamber, and three new leaves from each plant under different treatments were collected in liquid nitrogen for subsequent total RNA extraction.

Determination of Physiological and Biochemical Indexes
The chlorophyll extraction method involved the use of 95% ethanol as the solvent [9].A test tube containing 0.10 g of experimental material was sealed with a stopper, and then 10 mL of the 95% ethanol solution was added.The mixture was left to extract in darkness at room temperature for a duration of 24 h.Once the test material turned white, its color was compared using a spectrophotometer, with distilled water used as the reference point (zero).Simultaneously, the absorbance of the extracted solution at wavelengths of 665 nm, 649 nm, and 470 nm was measured to determine the chlorophyll a and b and carotenoid content.
The spectrophotometric assay of superoxide dismutase (SOD) activity was performed using the reduction of highly water-soluble tetrazolium salt (WST-8) by Xanthine-Xanthine oxidase [10].
The guaiacol colorimetry technique was employed to quantify peroxidase (POD) activity, utilizing absorbance changes at 470 nm as the basis for analysis [11].
The ultraviolet absorption method introduced by Johansson and Borg was used to determine catalase (CAT) activity [12].This method was based on the reaction of the enzyme with methanol in the presence of an optimal concentration of hydrogen peroxide.The formaldehyde produced was measured spectrophotometrically using 4-amino-3-hydrazino-5-mercapto-1,2,4-triazole (Purpald) as a chromogen.
The reduced glutathione (GSH) content was measured following Nishimoto et al. [13].Approximately 0.10 g of leaf sample was combined with 1 mL of extract solution, homogenized in an ice bath, and then centrifuged at 12,000 rpm at 4 • C for 10 min.The resulting supernatant was collected.Distilled water, the supernatant, and two reagents from a reduced glutathione kit were added to a 96-well plate.After thorough mixing and incubation for 2 min, the absorbance values at 412 nm of both the measurement tube and blank tube were measured to calculate the GSH content.
The amino acid content was assayed using the method of Chen et al. [14].Approximately 0.10 g of leaves was extracted with 10 mL of distilled water in a boiling water bath for 45 min.After centrifugation at 10,000× g for 10 min, 1 mL of supernatant was obtained and mixed with 0.5 mL of a 2% ninhydrin reagent and 0.5 mL of pH 8.0 phosphate buffer.The reaction mixture was heated in a boiling-water bath for 15 min.After cooling to room temperature, 8 mL of distilled water was added.The mixture stood for an additional 10 min before the absorbance was recorded at a wavelength of 562 nm using a MetashUV-5200 UV-vis spectrophotometer.
The bicinchoninic acid (BCA) assay, which involved two reactions, was used to determine the content of soluble protein [15].The first reaction resulted in the formation of a copper ion complex with amide bonds, leading to the reduction of copper in an alkaline solution.The second reaction involved the reduction of the BCA reagent primarily by the reduced copper-amide bond complex, as well as by tyrosine and tryptophan residues.

RNA Extraction and RNA-Seq Sequencing
The total RNA extraction was conducted using a plant total RNA extraction kit.The quality of the extracted RNA, including its degradation level and contamination status, were assessed through agarose gel electrophoresis.The concentration of RNA was accurately quantified using Qubit2.0(Thermo Scientific, Waltham, MA, USA), while the integrity of RNA was precisely evaluated with Agilent 2100 (Agilent Technologies Co., Ltd., Santa Clara, CA, USA).Furthermore, the purity of the obtained RNA was determined using Nanodrop spectrophotometer analysis.After meeting all qualification criteria for sample selection and preparation steps, library construction followed by Illumina sequencing were carried out [16].

Identification of Differential Gene Expression
In order to ensure the quality and reliability of data analysis, we used Trimmomatic software to filter the sequencing data.Specifically, we removed a subset of reads from the original data that contained adapters and exhibited low sequencing quality (Qphred ≤ 20 base number accounted for more than 50% of the entire read length) [17].Subsequently, Trinity was utilized to concatenate the clean reads, resulting in the generation of Unigene sequences, which served as reference sequences (Ref) [18].RSEM software was then applied to compare each sample's clean reads with Ref.Any reads with a comparative quality lower than 10 were filtered out accordingly.Additionally, any unpaired reads were matched to multiple regions of the genome.Finally, RSEM software was used to calculate the TPM value (transcripts per million kilobases) 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 the criterion q-value < 0.05 [16].

GO and KEGG Enrichment Analysis of DEGs
Enrichment analysis of Gene Ontology (GO) functional items in differentially expressed genes reveals significant correlations with biological functions compared to the genomic background.Firstly, all differentially expressed genes were mapped to each term in the GO database, and the gene count for each term was calculated [19].Subsequently, it was observed that there was significant enrichment of differentially expressed genes compared to the entire genome background.The enricher function of cluster-Profiler utilizes a hypergeometric distribution test as part of its GO enrichment analysis method to identify significantly enriched GO terms (Padj ≤ 0.05).
Pathway enrichment analysis enables 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 genome chemistry and system function data.The unit of analysis in significant pathway enrichment analysis was the KEGG pathway [20].A hypergeometric test was employed to identify pathways where differential genes were significantly enriched compared to all annotated genes.Similarly, a padj threshold less than 0.05 was used to determine significant enrichment in the KEGG pathway analysis.

Analysis of WGCNA
The weighted gene co-expression network analysis (WGCNA) algorithm was used to define the gene co-expression correlation matrix and the adjacency function of the gene network.Various coefficients were then computed for different nodes, and a hierarchical clustering tree was constructed to represent distinct branches of various gene modules [21].Hub genes were determined based on their connectivity degree within each module, followed by an exploration of the relationship between these hub genes and the effects of nano-selenium and organic selenium on leaves.

Quantitative Real-Time PCR
To validate the reliability of RNA-Seq results, we assessed nine differentially expressed genes (DEGs) identified with RNA-Seq using quantitative real-time polymerase chain reaction (qRT-PCR).We confirmed the expression patterns of these DEGs in leaves treated with nano-selenium and organic selenium.SiActin was used as the internal reference, and primers were designed based on the sequences of the DEGs.The RNA from the constructed library was reverse-transcribed into cDNA using PrimeScript TM 1st strand cDNA Synthesis Kit, followed by the incorporation of 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, which was then subjected to quantitative analysis [20].Finally, gene expression levels were calculated using the 2 −∆∆Ct method and analyzed [22].

Statistical Analysis
The data were sorted and graphs were created using Excel 2019 software.Statistical analysis was performed using SPSS 16.0 software (SPSS Inc., Chicago, IL, USA).Prior to statistical analysis, all data underwent normality and homogeneity tests.A one-way ANOVA was conducted to determine significant differences in physiological indexes of Suzhouqing under different organic selenium and nano-selenium supplementation treatments.The significance level for testing the differences was set at p < 0.05.

Photosynthetic Pigments
After spraying different concentrations of organic selenium and nano-selenium on Suzhouqing leaves, the chlorophyll content of the leaves was measured and recorded in Table 1.Compared to the control group, treatments with LO5, LO10, and LO20 showed an increase in chlorophyll a content by 16.85%, 6.74%, and 71.91%, respectively, with the LO20 treatment exhibiting the most significant effect.Similarly, LN5, LN10, and LN20 treatments resulted in an increase in chlorophyll a content by 17.98%, 15.73%, and 11.24%, respectively; among these treatments, LN10 had the highest chlorophyll a content.The application of LO5 led to a 4.44% increase in chlorophyll b content, while LO20 exhibited a significant increase of 40%.Furthermore, the carotenoid content significantly increased for all three treatments: LO5 (43.89%),LO10 (31.43%), and especially LO20 (1.31 times).However, there was no significant difference between the carotenoid contents at concentrations LN10 and LN20, indicating that effective increases could still be achieved even with lower concentrations.

Antioxidant Properties of Leaves
As illustrated in Table 2, the SOD activity of leaves treated with LO5, LO10, and LO20 was 1.73 times, 1.39 times, and 4.41 times that of the control group, respectively, with statistically significant differences observed; notably, the highest level was achieved under LO20 treatment conditions.In contrast, the SOD activity displayed a trend of initial decline followed by subsequent elevation within the LN5, LN10, and LN20 treatment groups.While the LN5 treatment resulted in a significant decrease of 39.87% in SOD activity compared to the control group, both the LN10 and LN20 treatments led to substantial increases of 31.39% and 38.80%, respectively; however, these differences were not statistically significant.The POD activity of Suzhouqing leaves treated with different concentrations of organic selenium is presented in Table 2.In the LO5, LO10, and LO20 treatment groups, the POD activity was 2.78 times, 1.27 times, and 2.00 times that of the control group, respectively.In conclusion, an appropriate concentration of organic selenium spray treatment could activate cellular antioxidant enzyme systems, enhance POD activity, improve stressresistance capabilities, and delay plant senescence.The table below illustrates the POD activity of leaves treated with different concentrations of nano-selenium.Compared to the control group, the leaves treated with LN5 exhibited a significant decrease in POD content by 20.01%.However, the leaves treated with LN10 and LN20 showed a significant increase in POD content of 24.00% and 31.69%,respectively; among them, LN20 reached its maximum value.
The LO20 treatment, as shown in Table 2, demonstrated the most significant enhancement of CAT activity, exhibiting an increase of 85.96% compared to the control group.The LO5 and LO10 treatments resulted in respective increases of 6.65% and 16.27% in CAT activity, with no significant difference observed for the LO5 treatment.These findings indicate that under appropriate concentrations of organic selenium, CAT can effectively inhibit lipid peroxidation and stimulate antioxidant responses.In the LN5, LN10, and LN20 treatment groups, CAT activity exhibited a pattern of initial decrease followed by an increase before decreasing again.Compared to the control group, leaves treated with LN5 showed a notable decrease in CAT content by 30.40%, reaching statistical significance.Conversely, leaves treated with LN10 and LN20 displayed significant increases in CAT content of 69.11% and 27.77%, respectively; among them, LN10 achieved the highest value.

Leaf Weight and Total Selenium Content
The leaf quality of individual plants is a crucial criterion for assessing the weight of leaves.As depicted in Table 3, in this experiment, the LO5, LO10, and LO20 treatments exhibited respective increases of 30.26%, 38.00%, and 0.95%.Similarly, the LN5, LN10, and LN20 treatments demonstrated increases of 16.19%, 28.66%, and 36.89%respectively.The total selenium content is presented in Table 3, showing that spraying organic selenium and nano-selenium led to an increase in the total selenium content.Specifically, the LO5, LO10, and LO20 treatments resulted in increases of 0.8 times, 1.9 times, and 1.4 times, respectively, while the LN5 treatment showed an increase of 90.00%, followed by LN10 with an increase of 2.20 times, and finally LN20 with an increase of 1.60 times.

Nutritional Quality of Leaves
Soluble amino acids are the fundamental constituents of an organism's protein.As depicted in Table 4, compared to the control group, the leaves treated with LO10 exhibited a significant decrease in amino acid content by 11.82%, while those treated with LO5 and LO20 showed significant increases of 19.09% and 32.27%, respectively, with LO20 reaching the highest value.LN5 treatment led to a significant increase in leaf amino acid content of 9.09%, whereas leaves treated with LN10 and LN20 experienced decreases of 7.73% and 4.55%, respectively, both reaching significance levels as well.The findings indicate that high concentrations of nano-selenium inhibit amino acid content.The soluble protein content in plants is often considered as a crucial indicator of plant oxidative senescence.A significant portion of the soluble protein found in mature plant leaves is ribulose diphosphate carboxylase (RuBPCase), which indirectly regulates RuBPCase levels, thereby enhancing photosynthetic capacity during the later stages of growth and effectively delaying crop senescence.As shown in Table 4, treatment with LO5, LO10, and LO20 resulted in increases of 17.51%, 22.46%, and 60.93%, respectively, in soluble protein content, with LO20 exhibiting the highest value at a statistically significant level.The increased soluble protein content also provided ample nitrogen sources for antioxidant enzyme synthesis, thereby enhancing the plant's ability to scavenge reactive oxygen species and delay leaf senescence while increasing the yield potential.In the LN5, LN10, and LN20 treatment groups, there was an initial increase followed by a subsequent decrease observed in the soluble protein levels over time.The leaves treated with LN5 showed an increase of 11.06% in soluble protein content, while those treated with LN10 exhibited a substantial increase of 49.39%.Additionally, the LN10 treatment demonstrated the highest level of soluble protein content among all treatments with significant differences.

Sequencing Results and Quality Analysis
The image data of the sequenced fragments were converted into sequence data (reads) using CASAVA base recognition.Transcriptomic libraries were obtained for the control group (LCK), nano-selenium treatment groups (LN5, LN10, and LN20), and organic selenium treatment groups (LO5, LO10, and LO20).The sequencing error rate for all samples was less than 0.1%, indicating high-quality results (Table 5).These findings demonstrate that RNA-Seq exhibited excellent quality with randomized and evenly distributed clean reads.By comparing reads that matched unique transcripts, differential genes could be identified.

Analysis of Differentially Expressed Genes
Correlation coefficients between groups were computed based on the TPM values of all genes in the sample.The results revealed that most Pearson correlation coefficients had an R 2 value greater than 0.80, indicating a strong correlation between gene expression levels across samples (Figure 1).

Sequencing Results and Quality Analysis
The image data of the sequenced fragments were converted into sequence data (reads) using CASAVA base recognition.Transcriptomic libraries were obtained for the control group (LCK), nano-selenium treatment groups (LN5, LN10, and LN20), and organic selenium treatment groups (LO5, LO10, and LO20).The sequencing error rate for all samples was less than 0.1%, indicating high-quality results (Table 5).These findings demonstrate that RNA-Seq exhibited excellent quality with randomized and evenly distributed clean reads.By comparing reads that matched unique transcripts, differential genes could be identified.Correlation coefficients between groups were computed based on the TPM values of all genes in the sample.The results revealed that most Pearson correlation coefficients had an R 2 value greater than 0.80, indicating a strong correlation between gene expression levels across samples (Figure 1).The volcano plot of differentially expressed genes (Figure 2) provides a visual representation of the overall distribution of these genes.The x-axis represents log2(Fold-Change), which indicates changes in gene expression multiples across different samples, The volcano plot of differentially expressed genes (Figure 2) provides a visual representation of the overall distribution of these genes.The x-axis represents log 2 (FoldChange), which indicates changes in gene expression multiples across different samples, while the y-axis represents log 10 (padj), which reflects the significance levels of expression differences.

GO Annotation and Enrichment Analysis of Differentially Expressed Genes
GO consists of three categories: biological process, cellular component, and molecular function.Figure 3 shows the number of up-regulated genes, down-regulated genes, and total genes within each GO category.In Figure 3a, metabolic process, plastid nucleoid, and alcohol binding had the highest number of up-regulated genes.In Figure 3b, chlorophyll metabolic process and regulation of dephosphorylation had a higher number of up-regulated genes along with tubulin complex and alcohol binding.In Figure 3c, regulation of dephosphorylation, cajal body, and tubulin complex along with alcohol binding showed a higher count for up-regulated genes.Similarly, in Figure 3d, metabolic process along with tubulin complex and alcohol binding displayed a higher count for up-regulated genes.In Figure 3e, metabolic process together with plastid nucleoid and phosphatase regulator activity exhibited a higher number of up-regulated genes.Lastly, in Figure 3f, sesquiterpenoid metabolic process demonstrated the largest count of up-regulated genes followed by protein phosphatase type 2A complex and protein histidine kinase binding.
Agriculture 2023, 13, x FOR PEER REVIEW 9 of 27 while the y-axis represents log10(padj), which reflects the significance levels of expression differences.

GO Annotation and Enrichment Analysis of Differentially Expressed Genes
GO consists of three categories: biological process, cellular component, and molecular function.Figure 3 shows the number of up-regulated genes, down-regulated genes, and total genes within each GO category.In Figure 3a, metabolic process, plastid nucleoid,    The pathway map illustrating the enrichment of differentially expressed genes is presented in Figure 5.To visualize the distribution of these genes within the pathway map, differential genes were annotated, with red indicating up-regulated genes and green representing down-regulated genes.
The light-harvesting protein complex in the pathway of light-antenna proteins comprises both up-regulated and down-regulated genes.The chlorophyll-binding subunits of photosystems I and II serve as internal antenna light-harvesting proteins for oxygenic photosynthesis.Phycobilisomes contain antenna proteins, while green plants possess lightharvesting chlorophyll protein complexes acting as peripheral antenna systems to enhance the efficiency of light energy absorption.
The expression of LHca2 (light-harvesting complex I chlorophyll a/b binding protein 2) was up-regulated in LN5 compared to LCK.Additionally, LHcb1 (light-harvesting complex II chlorophyll a/b binding protein 1) showed an increase in expression, as well as LHcb5 (light-harvesting complex II chlorophyll a/b binding protein 5).In LN10 vs. LCK, there was an up-regulation of LHca1 (light-harvesting complex I chlorophyll a/b binding protein 1), along with the up-regulation of LHca2 and LHcb7 (light-harvesting complex II chlorophyll a/b binding protein 7).Similarly, in LN20 vs. LCK, there was an up-regulation of LHca1, along with the up-regulation of LHca2 and LHcb4 (light-trapping complex II chlorophyll a/b binding protein).Furthermore, in LO5 vs. LCK, there was an up-regulation of LHca1.Specifically, there was an increase in the expression of chlorophyll a/b binding protein 2, also known as LHca2.Moreover, both LHcb1 and Lhcb2 were found to be up-regulated.Lastly, both LHcb4 and LHcb6 exhibited increased expression levels.LO10 vs. LCK specifically showed an increase in the expression level of LHca1 while also showing increased levels of LHcb5.

Weighted Gene Co-Expression Network Analysis (WGCNA)
In order to identify gene co-expression modules, WGCNA was conducted for all expressed genes.Firstly, the co-expression network was constructed based on the transcriptome sequencing results using a defined fitting index of R 2 = 0.80.Low-expression genes were excluded and the cluster of co-expressed genes was dynamically cut.A cluster tree was then generated according to the correlation between gene expression levels (Figure 6), resulting in module divisions.Each color in the figure represents a specific module that includes genes belonging to the same cluster in the clustering tree.If certain genes consistently exhibit similar changes in expression during a physiological process or across different tissues, it suggests their functional relationship and they can be defined as a module.The vertical distance reflects the dissimilarity between two nodes (genes).The correlation heat map between modules (Figure 7) can be divided into two parts: the upper part clusters modules based on their eigengene values, while each row and column in the bottom half represents an individual module.Blue indicates a negative correlation with a coefficient ranging from 0 to 0.5; values closer to zero indicate stronger negative correlations.Red indicates a positive correlation with coefficients ranging from 0.5 to 1; values closer to one indicate stronger positive correlations.The results demonstrated a significant correlation of 0.75 between the MECoral module and MEDarkturquoise module, thereby establishing these two modules as potential target gene modules.Subsequently, an independent analysis was conducted on the expression levels of all genes and eigenvector genes within the two modules across all samples.The gene expression levels within each module showed a significant degree of correlation.And the expression levels of eigenvector genes were highly related to the overall expression levels of their respective modules.This finding suggests that eigenvector genes in target modules can effectively represent the overall gene expressions within their corresponding modules (Figure 8).Finally, hub genes from the two modules were selected based on their connectivity ranking with other genes; specifically, those with higher kWithin rankings indicating central hub positions were chosen.The MEDarkturquoise module yielded five core genes, while the MECoral module identified four core genes, resulting in a total of nine hub genes (LOC103857346, Actin-related protein 2/3 complex subunit 1A; LOC103828797, Actin-related protein 2/3 complex subunit 2A; LOC103874199, Superoxide dismutase [Fe] 3; LOC103838242, Monodehydroascorbate reductase; LOC103846588, Thioredoxin M2; LOC103873206, ABC transporter A family member 7; LOC106353797, ABC transporter B family member 26; LOC103859979, ATP sulfurylase 1; LOC103859070, Cysteine synthase D1) associated with the antioxidant stress, absorption, transport, and metabolism of selenium.All the above hub genes were related to endocytosis, antioxidant activity, the ABC transporter, and selenium metabolism.Further analysis showed that the functional enrichment analysis results of the nine hub genes were similar to the results of the GO enrichment analysis and KEGG pathway analysis of DEGs mentioned in parts 3.5.3and 3.5.4,respectively, indicating that the above analysis was reliable.Real-time fluorescence quantitative PCR was conducted to validate these nine hub genes, as depicted in Figure 9.
Agriculture 2023, 13, x FOR PEER REVIEW 10 of 27 of dephosphorylation, cajal body, and tubulin complex along with alcohol binding showed a higher count for up-regulated genes.Similarly, in Figure 3d, metabolic process along with tubulin complex and alcohol binding displayed a higher count for up-regulated genes.In Figure 3e, metabolic process together with plastid nucleoid and phosphatase regulator activity exhibited a higher number of up-regulated genes.Lastly, in Figure 3f, sesquiterpenoid metabolic process demonstrated the largest count of up-regulated genes followed by protein phosphatase type 2A complex and protein histidine kinase binding.(

2) Enrichment of KEGG pathway diagram
The pathway map illustrating the enrichment of differentially expressed genes is presented in Figure 5.To visualize the distribution of these genes within the pathway map, differential genes were annotated, with red indicating up-regulated genes and green representing down-regulated genes.
The light-harvesting protein complex in the pathway of light-antenna proteins comprises both up-regulated and down-regulated genes.The chlorophyll-binding subunits of photosystems I and II serve as internal antenna light-harvesting proteins for oxygenic photosynthesis.Phycobilisomes contain antenna proteins, while green plants possess light-harvesting chlorophyll protein complexes acting as peripheral antenna systems to enhance the efficiency of light energy absorption.
The expression of LHca2 (light-harvesting complex I chlorophyll a/b binding protein 2) was up-regulated in LN5 compared to LCK.Additionally, LHcb1 (light-harvesting complex II chlorophyll a/b binding protein 1) showed an increase in expression, as well as LHcb5 (light-harvesting complex II chlorophyll a/b binding protein 5).In LN10 vs. LCK, there was an up-regulation of LHca1 (light-harvesting complex I chlorophyll a/b binding protein 1), along with the up-regulation of LHca2 and LHcb7 (light-harvesting complex II chlorophyll a/b binding protein 7).Similarly, in LN20 vs. LCK, there was an up-regulation of LHca1, along with the up-regulation of LHca2 and LHcb4 (light-trapping complex II chlorophyll a/b binding protein).Furthermore, in LO5 vs. LCK, there was an up-regulation of LHca1.Specifically, there was an increase in the expression of chlorophyll a/b binding protein 2, also known as LHca2.Moreover, both LHcb1 and Lhcb2 were found to be up-regulated.Lastly, both LHcb4 and LHcb6 exhibited increased expression levels.LO10 vs. LCK specifically showed an increase in the expression level of LHca1 while also showing increased levels of LHcb5.In order to identify gene co-expression modules, WGCNA was conducted for all expressed genes.Firstly, the co-expression network was constructed based on the transcriptome sequencing results using a defined fitting index of R 2 = 0.80.Low-expression genes were excluded and the cluster of co-expressed genes was dynamically cut.A cluster tree was then generated according to the correlation between gene expression levels (Figure hub genes were related to endocytosis, antioxidant activity, the ABC transporter, and selenium metabolism.Further analysis showed that the functional enrichment analysis results of the nine hub genes were similar to the results of the GO enrichment analysis and KEGG pathway analysis of DEGs mentioned in parts 3.5.3and 3.5.4,respectively, indicating that the above analysis was reliable.Real-time fluorescence quantitative PCR was conducted to validate these nine hub genes, as depicted in Figure 9.

Fluorescence Quantitative PCR Analysis
The differential genes were validated using qRT-PCR, and the mapping analysis yielded a high R 2 value calculated based on the 2 −∆∆Ct method.The concordance between the qRT-PCR validation and RNA-Seq results demonstrates the reliability of the experimental findings (Figure 9).3.5.7.Differential Gene Function Analysis (1) Genes related to selenium absorption-endocytosis Endocytosis is a cellular mechanism employed for the acquisition of extracellular nutrients.Further analysis showed that endocytosis included not only the hub genes

Fluorescence Quantitative PCR Analysis
The differential genes were validated using qRT-PCR, and the mapping analysis yielded a high R 2 value calculated based on the 2 −∆∆Ct method.The concordance between the qRT-PCR validation and RNA-Seq results demonstrates the reliability of the experimental findings (Figure 9).

Differential Gene Function Analysis
(1) Genes related to selenium absorption-endocytosis Endocytosis is a cellular mechanism employed for the acquisition of extracellular nutrients.Further analysis showed that endocytosis included not only the hub genes (LOC103857346, Actin-related protein 2/3 complex subunit 1A; LOC103828797, Actinrelated protein 2/3 complex subunit 2A), but also other genes, as shown in Table 6.Transcriptome analysis revealed 32 differentially expressed genes associated with endocytic proteins, as presented in Table 6.Notably, Actin-related protein 2/3 complex subunit 1A, nuclear-pore anchor, and vacuolar protein sorting-associated protein exhibited enhanced expression upon treatment with organic selenium and nano-selenium 2 homolog 3, along with up-regulation observed in vacuolar protein sorting-associated protein 28 homolog 2 and Vacuolar protein sorting-associated protein 60.This suggests that endocytosis may facilitate the absorption of selenium.(2) Genes involved in antioxidant stress Low temperature can induce the generation of reactive oxygen species (ROS) in plants, and excessive ROS accumulation can trigger cell apoptosis.To enhance the ROS response, plants up-regulate the expression of antioxidant coding genes, promote antioxidants biosynthesis, and thereby improve their antioxidant capacity.In addition to the hub genes mentioned above, there were other genes involved in antioxidant activity.In both the organic selenium and nano-selenium treatment groups, superoxide dismutase (SOD), superoxide dismutase [Fe] 3, superoxide dismutase [Fe] 2, monodehydroascorbate reductase, and thioredoxin M2 were significantly up-regulated, playing a crucial role in combating oxidative stress caused by low temperature in leaves (Table 7).(3) Selenium transport-related gene-ABC transporter The ABC transporter plays a crucial role in the cell membrane by facilitating the transport of substances across the cell or extracellular space through ATP.As shown in Table 8, the comparative analysis revealed differential expression of 23 ABC transporter genes, with up-regulation observed for ABC transporter B family member 29, ABC transporter C family member 10, ABC transporter E family member 2, ABC transporter G family member 11, ABC transporter G family member 40, and ABC transporter I family member 1. (4) Genes related to selenium metabolism The metabolic pathways of cysteine and methionine, sulfur, selenide, and glutathione, which are closely associated with selenium metabolism, were analyzed.The subsequent analysis revealed that selenium metabolism involved not only the hub genes but also other genes, as indicated in Table 9.After selenium absorption, the expression of genes related to the sulfur metabolic pathway was either up-regulated or down-regulated.The expression of chloroplastic ATP-ylase 1 genes was up-regulated following treatment with organic selenium and nano-selenium.Hub genes controlling the expression of cysteine synthase showed varying degrees of up-regulation.Furthermore, it was observed that genes involved in selenide metabolism-related enzymes were up-regulated.These included glutathione S-transferase U12 and S-sulfo-L-cysteine synthase (O-acetyl-L-serine-dependent).On the other hand, glutathionyl-hydroquinone reductase YqjG and probable glutathione peroxidase 8 were down-regulated.In summary, there was a higher overall abundance of up-regulated gene expression compared to down-regulated gene expression, and significant differences in gene expression levels were observed between organic selenium and nano-selenium treatments.Notably, under nano-selenium treatment conditions, the hub genes exhibited significantly higher expression levels than under organic selenium treatment conditions, indicating a superior absorption effect for nano-selenium.Low temperature stress can cause damage to chloroplasts in plants, which directly affects the production of photosynthetic pigments and the process of photosynthesis [23].Zhang et al. used three different strains of Dendrobium candidum seedlings as materials and sprayed varying concentrations of sodium selenite solution on the leaves at low temperatures [24].The results showed that the proper selenium treatment concentration (0.10-0.20 mg L −1 ) could alleviate the decrease in chlorophyll content.It was speculated that selenium may promote the absorption of mineral elements related to chlorophyll synthesis in plants or participate in the synthesis of chlorophyll precursors in the form of seleno-amino acids.Sun studied the mitigating effect of exogenous selenium on strawberry seedlings under low temperature stress and its influence on the ASA-GSH cycle.They found that after spraying exogenous selenium, it effectively alleviated chlorophyll degradation caused by low temperature stress in strawberry seedlings and significantly increased the chlorophyll a/b value.This might be attributed to selenium's ability to effectively reduce the accumulation levels of protoporphyrin-IX and MG-protoporphyrin esters related to chlorophyll synthesis in plants [25].
Chen found that the application of different valences of selenium could significantly increase the content of photosynthetic pigments in tobacco leaves under low temperature stress.The content of chlorophyll a and chlorophyll b was higher under sodium selenate treatment compared to sodium selenite treatment [26].Our results are consistent with the above findings.Under low temperature stress conditions, the application of appropriate concentrations of organic selenium (20 mg L −1 ) and nano-selenium (10 mg L −1 ) significantly increased the chlorophyll content and chlorophyll a/b ratio in leaves.This is because it up-regulates the expression of light-harvesting complex I chlorophyll a/b binding proteins, thereby enhancing plant absorption, transfer, and utilization of photons to improve light energy conversion efficiency.

Effects of Spraying Organic Selenium and Nano-Selenium Fertilizer on the Regulation of Antioxidant System of Leaves under Low Temperature Stress
Our findings reveal that plants utilize a collaborative and mutually influential mechanism involving multiple enzymes to protect cellular structure and function from damage caused by reactive oxygen species accumulation under low temperature stress.SOD serves as the primary defense against superoxide anions by catalyzing their dismutation into oxygen and hydrogen peroxide.Additionally, CAT and POD further decompose H 2 O 2 to prevent oxidative damage to plant cells.The activity of SOD in leaves was higher than that of CAT and POD in this study.It can be speculated that when the plant encounters low temperatures, SOD reacts first, converting the superoxide anion into H 2 O 2 .Moreover, exogenous application of selenium significantly enhances SOD activity by up-regulating the transcriptional expression of SOD genes, including superoxide dismutase [Fe] 2 and superoxide dismutase [Fe] 3. The increased SOD activity leads to elevated H 2 O 2 production.However, selenium also notably improves the activities of CAT and POD enzymes for effectively scavenging the generated H 2 O 2 .Song et al. demonstrated that lettuce seedlings could mitigate membrane lipid peroxidation damage under low temperature stress by increasing the activities of SOD and POD [27].Similarly, Wang et al. revealed that suitable spraying concentrations of nano-selenium fertilizer enhanced SOD and POD activities [28].We obtained similar results to those in the literature mentioned above.Additionally, our findings were consistent with previous studies conducted on strawberry, tobacco, and wheat that reported an increase in POD activity due to selenium supplementation [25,26,29].However, Liu et al. observed a certain inhibitory effect of selenium on POD activity during their investigation on Pyrus bretschneider Rehd.'Dangshan Suli' [30].This discrepancy may be attributed to selenium's ability to enhance antioxidant enzyme activities within SOD, CAT, and AsA-GSH while reducing substrate availability for POD enzymatic reactions.

Effects of Spraying Organic Selenium and Nano-Selenium Fertilizer on Leaf Nutritional Quality under Low Temperature Stress
Amino acids in plants can regulate plant growth, promote nutrient absorption, and enhance metabolic function.Liu demonstrated that under low temperature stress, the content of free amino acids in the exogenous selenium treatment group was 16.30% higher than that in the control group.This indicated that the application of exogenous selenium effectively increased the content of free amino acids in tea trees [31].Cao found that under low-temperature conditions, the soluble protein content of potato leaves exhibited an increasing trend with an increase in foliar spray concentration [32].The soluble protein content of potato leaves with a foliar spray selenium concentration of 20 mg L −1 was 20.90 mg L −1 , which was 19.81% higher than that of the control group.These results indicated that spraying selenium on leaf surfaces could increase the soluble protein content to some extent and improve cold tolerance in potato plants.The findings from this study were consistent with those mentioned above.Nano-selenium spraying could elevate the soluble protein content by directly involving selenocysteine (Se-Cys) and selenomethionine (Se-Met) in protein synthesis while reducing cysteine (Cys) and methionine (Met) levels among free amino acids.

Absorption of Organic Selenium and Nano-Selenium at Different Concentrations in Leaves
Currently, research is primarily focused on the mechanism of inorganic selenium absorption in plant roots, with limited reports on organic selenium and nano-selenium absorption in leaves.It has been hypothesized that transporters related to plant cysteine and methionine may be involved in some processes of organic selenium absorption.Kong et al. discovered significant differences in genes related to the endocytosis process after treatment with nano-selenium compared to the control group [33].Therefore, it was hypothesized that millet's absorption of nano-selenium was linked to the endocytosis process.And further tests were needed to verify whether genes related to PIP5K, Hsc70, PLD, and ArfGAP were involved in millet's selenium uptake and how they participate in this process.This study concluded that absorption of nano-selenium and organic selenium was associated with endocytosis, aligning with previous findings.

Transport of Different Concentrations of Organic Selenium and Nano-Selenium in Leaves
This study revealed that ABC transporters play a crucial role in promoting plant growth by facilitating the transportation of nutrients, hormones, lipids, secondary metabolites, metal ions, and foreign substances.Additionally, they contribute to disease resistance and response to stress in plants.Moreover, these transporters are involved in regulating ion channels.The ABCB-subgroup transporters primarily function in the transportation of heavy metals and secondary metabolites for resistance purposes.On the other hand, the ABCC subfamily is responsible for transporting glutathione compounds.Byrne et al. reported that 28 ABC transporters in ryegrass might be associated with selenide accumulation [34].In the investigation of genes related to selenium absorption and metabolism in tea tree roots conducted by Hu, 16 hub genes of ABC transporter were identified, suggesting their potential involvement in selenate transport within roots [35].However, Kong et al. examined selenium-response genes in millet leaves and found that out of seven differential genes belonging to ABC transporter family, five were up-regulated, while two were down-regulated, contradicting previous studies' findings [26].These discrepancies might arise from variations observed among different species regarding the expression patterns of ABC transporters.

Metabolism of Different Concentrations of Organic Selenium and Nano-Selenium in Leaves
Usually, selenite is reduced to selenide by sulfite reductase or glutathione, and selenide can also be produced through various mechanisms such as methionol in the action of selenium-binding proteins.Subsequently, cysteine synthase catalyzes the conversion of selenide into selenocysteine (SeCys) [36].Kong et al. discovered that hub genes regulating the expression of Se-binding protein 1 and cysteine synthase were up-regulated to varying degrees during this process [33].The hub genes associated with cysteine synthase in this study exhibit similarities to those identified in previous research.Following the biosynthesis of selenocysteine, a portion of it undergoes further conversion into selenocysta by cystathionine-gamma-synthetase (CγS), which was subsequently decomposed into se-homocysteine (SeHCys) by cystathionine-β-lyase (CBL).Finally, under the influence of cysteine methyltransferase (MET), it was transformed into selenomethionine (SeMet).Some synthetic SeCys and SeMet could substitute Cys and Met within plant protein chains, thereby influencing their structure and function [37].Other non-protein selenium amino acids like methylselenol and its derivatives, as well as non-amino acid selenides, could be generated through the action of methionin-gamma-lyase to mitigate the toxic effects exerted by SeMet on plants' growth processes [38,39].Inserting the SMT gene in Arabidopsis thaliana resulted in an increased content of SeMSeCys along with enhanced tolerance towards selenium [40].Hu did not detect any SMT genes in tea tree roots following selenium treatment, indicating that SMT was not implicated in selenium regulation [35].Similarly, Tan found no differential expression of SMT due to selenium induction [41].In this study, while no SMT was identified, the crucial gene homocysteine S-methyltransferase was discovered, and it was observed that the selenocysteine S-methyltransferase gene CsSMT mRNA could interact with the homocysteine S-methyltransferase gene of Arabidopsis thaliana [42].Furthermore, the up-regulation of glutaredoxin-C1, glutaredoxin-C6, and glutathione S-transferase U12 suggested that selenium played a role in regulating the synthesis of glutathione, cysteine, and methionine, which ultimately affects its accumulation.

Conclusions
The impact of organic selenium and nano-selenium at various concentrations on Suzhouqing's growth, nutrition, and antioxidant properties under low temperatures was examined.The findings demonstrated that exogenous selenium application significantly enhanced leaf photosynthesis, enhanced antioxidant activity, increased leaf weight, and promoted selenium uptake and translocation within leaves.The transcriptome sequencing of leaves under different treatments was conducted using RNA-Seq technology, and the specific reasons were systematically analyzed.The study revealed that selenium influenced the leaf weight and total selenium content through three main mechanisms.First, it protected photosynthetic pigments and boosted photosynthetic capacity by up-regulating LHca2, LHcb1, LHca1, and LHcb4.Second, it enhanced antioxidant capacity by elevating the expression of genes such as superoxide dismutase and monodehydroascorbic acid.Third, it facilitated selenium absorption through endocytosis, transported selenium using the ABC transport gene family, and regulated selenium metabolism-related genes like cysteine synthetase and glutaredoxin.Subsequently, nine hub genes identified with WGCNA were closely associated with the mentioned processes of endocytosis, antioxidant activity, ABC

Figure 1 .
Figure 1.. Heat map of correlation coefficients between samples.

Figure 1 .
Figure 1.Heat map of correlation coefficients between samples.

Figure 2 .
Figure 2. Volcano map of differentially expressed genes ("−1" indicates down-regulated, "0" indicates unchanged and "1" indicates up-regulated).Note: The horizontal coordinate is log2(fold change), that is, the pair value of the fold change value; The ordinate is -log10 (padj), which is the inverse of the logarithm of the corrected p-value.The up-regulated genes are indicated by red dots, while the down-regulated genes are represented by blue dots.

Figure 2 .
Figure 2. Volcano map of differentially expressed genes ("−1" indicates down-regulated, "0" indicates unchanged and "1" indicates up-regulated).Note: The horizontal coordinate is log 2 (fold change), that is, the pair value of the fold change value; The ordinate is -log 10 (padj), which is the inverse of the logarithm of the corrected p-value.The up-regulated genes are indicated by red dots, while the down-regulated genes are represented by blue dots.3.5.4.KEGG Annotation and Enrichment Analysis of Differentially Expressed Genes (1) KEGG annotation of differentially expressed genes Figure 4 illustrates the proportion of annotated genes and included pathways by comparing Unigene with the KEGG database.

Figure 4 .
Figure 4. KEGG metabolic pathway classification statistics.Note: The ordinate is the annotated KEGG pathway, and the abscissa is the proportion of genes annotated to the pathway.

Figure 4 . 27 Figure 5 .
Figure 4. KEGG metabolic pathway classification statistics.Note: The ordinate is the annotated KEGG pathway, and the abscissa is the proportion of genes annotated to the pathway.Agriculture 2023, 13, x FOR PEER REVIEW 14 of 27

Figure 6 .
Figure 6.Hierarchical clusterting tree.Note:The color coding indicates that each gene in the clustering tree belongs to a specific module.A module is defined as a group of genes with similar expression changes during a physiological process.In the upper part of the tree, the vertical distance represents the genetic distance between two nodes.

Figure 6 .
Figure 6.Hierarchical clusterting tree.Note:The color coding indicates that each gene in the clustering tree belongs to a specific module.A module is defined as a group of genes with similar expression changes during a physiological process.In the upper part of the tree, the vertical distance represents the genetic distance between two nodes.

Figure 7 .
Figure 7. Module eigenegene correlation between different modules of WGCNA.Note: The correlation heat map between modules can be divided into two parts, with the upper part clustering modules according to their eigengene values.The ordinate represents the degree of dissimilarity of nodes; in the bottom half of the graph, each row and column represents a module.The darker (redder) the square color, the stronger the correlation; the lighter the square color, the weaker the correlation.

Figure 7 .
Figure 7. Module eigenegene correlation between different modules of WGCNA.Note: The correlation heat map between modules can be divided into two parts, with the upper part clustering modules according to their eigengene values.The ordinate represents the degree of dissimilarity of nodes; in the bottom half of the graph, each row and column represents a module.The darker (redder) the square color, the stronger the correlation; the lighter the square color, the weaker the correlation.

Figure 8 .
Figure 8. Expression levels of all genes and corresponding module eigengene in MEdarkturquoise (left) module and MECoral (right) module.Note: The module gene expression pattern map can be divided into two parts: the header annotation is the module name, the left side of the above figure is the gene name, the horizontal coordinate of the above figure is the sample name, and the above figure is the expression calorimetric map of the gene in the module in different samples, with red indicating high expression and green indicating low expression.The following figure shows the expression patterns of module eigenvalues in different samples, and the abscissa is the sample name.

Figure 8 .
Figure 8. Expression levels of all genes and corresponding module eigengene in MEdarkturquoise (left) module and MECoral (right) module.Note: The module gene expression pattern map can be divided into two parts: the header annotation is the module name, the left side of the above figure is the gene name, the horizontal coordinate of the above figure is the sample name, and the above figure is the expression calorimetric map of the gene in the module in different samples, with red indicating high expression and green indicating low expression.The following figure shows the expression patterns of module eigenvalues in different samples, and the abscissa is the sample name.Agriculture 2023, 13, x FOR PEER REVIEW 18 of 27

Table 1 .
Photosynthetic pigment content of Suzhouqing treated with different concentrations of organic

Table 2 .
Antioxidant properties of green Suzhouqing treated with different concentrations of organic selenium and nano-selenium.

Table 3 .
Weight and total selenium content of Suzhouqing leaves.
Note: Lowercase letters indicate significant differences between treatments (p < 0.05) in this paper.

Table 4 .
Nutritional quality of Suzhouqing treated with different concentrations of organic selenium and nano-selenium.
Note: Lowercase letters indicate significant differences between treatments (p < 0.05) in this paper.

Table 6 .
The related genes in endocytosis.

Table 7 .
Genes involved in antioxidative stress.

Table 8 .
The genes of the ABC transporter gene family.