The 9930 cultivated cucumber inbred lines, typical of North China fresh market type cucumber, were used for this study. Cucumbers were grown in growth chambers under 16 h light/8 h dark at 28°C day/18°C night to S fumigation and control treatment, respectively. There were 20 cucumber seedlings in a growth chamber. The cucumber plants were S fumigated using 2 h/day, 0.02 g/day/m³ at 40-days-old, where cucumber seedlings entered the rapid growth phase (Gómez et al., 2003).

Physiological, transcriptomic, and metabolomic data were collected from three biological duplicates in cucumber fruit and leaf with developmental stages. Each biologically duplicated sample pool was collected from five cucumber plants. The treatments included ovary (0 days post-pollination, F0), commercial fruit (12 days post-pollination, F12), young leaf (ovary derived leaf, L0), and mature leaf (commercial fruit-derived leaf, L12). Supplementary Table 1 shows the specific information of sampling sites.

The chlorophyll (Chl) contents and Chl fluorescence parameters were measured in non-fumigated and S-fumigated L0 and L12, respectively. In brief, the Chl were extracted from 1.0 g of the fresh leaves using 95% ethanol for 24 h. The Chl extracts of at least three replicates were measured at 470, 649, and 665 nm wavelengths using the SPECORD 210 PLUS spectrophotometer (Analytik Jena GmbH, Jena, Germany).

Leaves were dark-adapted for 30 min using a pulse amplitude modulated Mini-PAM fluorometer (Heinz Walz, Effeltrich, Germany), and Chl fluorescence parameters were measured in eight replicates as previously described (Li et al., 2021). The minimum (F0) and maximum fluorescence (Fm) of dark-adapted leaves, actual (F), and maximum fluorescence of light-adapted leaves (Fm′) were evaluated. The variable to maximum fluorescence ratio (Fv/Fm) was calculated as (Fm − F0)/Fm, whereas the effective quantum yield (ΦPSII) was calculated as (Fm′ − F)/Fm′.

The activities of antioxidant enzymes, i.e., superoxide dismutase (SOD) and peroxidase (POD), and contents of the organic compound malondialdehyde (MDA) were determined in L0 and L12 under non-fumigated and S-fumigated treatments, respectively. The activities of SOD and POD were determined using a superoxide dismutase assay kit (Solarbio, Beijing, China) and a peroxidase assay kit (Solarbio, Beijing, China), respectively. The malondialdehyde assay kit (Solarbio, Beijing, China) measured the MDA contents. In brief, collected leaves were ground using liquid nitrogen, and SOD, POD, and MDA were extracted using the extraction buffer. The extracts were centrifuged at 12,000 × g for 10 min, and the supernatants were then measured for SOD, POD activities, and MDA contents using the multimode microplate reader (Tecan, Switzerland) (Yin et al., 2018). Each indicator had three replicate measurements.

Fruits were collected from 0, 3, 6, 9, and 12 days post-pollination to measure fruit development using the ratio of fruit length (L) and diameter (D). The fruit quality was evaluated by the fruit mass, peel color, and total soluble solids (TSS, %). At F12, fruit mass was measured with 30 duplicates using an analytical balance (0.01 g accuracy). The peel color of eight fruits was determined using a CR-10 Plus model colorimeter (Minolta Camera, Osaka, Japan) in the CIELAB color space. Meanwhile, the TSS (%) was determined by three fruits using a digital handheld PAL-1 pocket refractometer (Atago, Tokyo, Japan) (Yun et al., 2013).

The total RNA was extracted from three biological replicates of L0, L12, F0, and F12 using the mirVana miRNA Isolation Kit (Ambion, Austin, TX, United States). The RNA integrity was determined using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). Then, cDNA libraries were prepared from 4 μg of total RNA using the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, United States). The prepared cDNA libraries were sequenced on the Illumina HiSeq X Ten platform, and 150 bp paired-end reads were generated by the OE Biotech Co., Ltd. (Shanghai, China). Sequenced reads were quality assessed, and low-quality sequences were removed. Cleaned reads were mapped to the Chinese long cucumber genome version 3¹ using HISAT2 (Kim et al., 2015). The read counts were estimated as the number of fragments per kilobase of the transcript sequences per million base pairs mapped (FPKM) values of expressed genes. The differentially expressed genes (DEGs) were identified using DESeq (2012) in R package (Anders and Huber, 2013) with false discovery rate (FDR) < 0.05 and | log2(Fold Change)| > = 1 (Ni et al., 2021).

Gene ontology (GO) enrichment analysis of DEGs was performed using R (version 4.0.4) based on the hypergeometric distribution. The weighted gene co-expression network analysis (WGCNA) was performed using R with default parameters to identify co-expressed gene modules. The adjacency matrix was constructed using Pearson’s coefficients between genes and transformed into the topological overlap matrix. Clustered genes were grouped into modules according to the expression patterns. The WGCNA package associated the phenotype data with the classified modules. The gene significance (GS) and module membership (MM) of each node and the weight of node-to-node were calculated using the WGCNA package in each module. The genes of the top 10% GS and MM combined with differential expression were selected. The information on these genes was extracted from the cucumber (Chinese Long) version 3 genome. Then, hub genes within each network were obtained according to the maximal clique centrality (MCC) method, and the module network data were visualized using Cytoscape version 3.7.1 (Shanno et al., 2003).

The gene expression of three independent biological replicates per sample was determined using quantitative real-time PCR (RT-qPCR). The total RNA was extracted using the RNAprep Pure Plant Kit (TIANGEN, Beijing, China), and cDNA was synthesized using a PrimeScript^TM RT Reagent Kit (Takara, Tokyo, Japan) with reverse transcriptase, following the instructions of the manufacturer. The hub genes were selected for RT-qPCR validation as previously described (Qin et al., 2021). The reactions were normalized using the cucumber actin gene (CuGenDB name: CsaV3_6G041900). Specific primer pairs were designed for specific unigenes and actin (Supplementary Table 2). The RNA level was expressed relative to the actin gene expression level following the 2^–ΔΔCT method.

The same RNA-seq samples from L0, L12, F0, and F12 were used for metabolite analysis. The OE Biotech Co., Ltd. (Shanghai, China) prepared samples, extracted, identified, and quantified metabolites. In brief, ∼80 mg of tissue was transferred into an Eppendorf tube containing 20 μl internal standard (2-chloro-L-phenylalanine in methanol, 0.3 mg/ml) and 1 ml extraction solvent with methanol/H₂O (4/1, v/v). The mixture was ground using an automatic high-throughput tissue crusher at 60 Hz, −20°C for 2 min. The extract was sonicated at 50 Hz, −20°C for 30 min, centrifuged at 13,000 rpm, and 4°C for 10 min.

The supernatant was analyzed using the ultra-performance liquid chromatography (UPLC) Triple TOF system (AB SCIEX, Framingham, MA, United States). The positive and negative modes were the ACQUITY UPLC BEH C18 column (1.7 μm, 2.1 × 100 mm) (Waters Corporation, Milford, MA, United States). Water (containing 0.1% formic acid) was the mobile phase A, and acetonitrile (containing 0.1% formic acid, v/v) was the mobile phase B. The flow rate was 0.4 ml/min, and the column temperature was 45°C. The quality control (QC) samples were used for internal standard (every 10 samples) throughout the analytical run. The ultra-performance liquid chromatography and tandem mass spectrometry (UPLC-MS/MS) raw data were analyzed using the Progenesis QI software (Waters Corporation, Milford, MA, United States), based on public^{2 3} and self-built databases.

The metabolic alterations among experimental groups were visualized using the principal component analysis (PCA) and (orthogonal) partial least-squares-discriminant analysis [(O)PLS-DA] via online resources⁴. The ellipse was defined as 95% CI of the modeled variation in the score plot. The variable importance in the projection (VIP) ranked the overall contribution of each variable to the OPLS-DA model. Variables with VIP > 1.0 were considered relevant for group discrimination. The differential metabolites were selected based on the combination of statistically significant threshold VIP values from the OPLS-DA model and the P-values from the two-tailed Student’s t-test of normalized peak areas. Metabolites with VIP > 1.0 and P < 0.05 were considered significantly changed. Meanwhile, the bidirectional orthogonal projections to latent structures (O2PLS) were used to integrate metabolome and transcriptome data via online resources⁵.

There were no obvious phenotypic changes in cucumber seedlings at 40 days after sulfur fumigation (Figure 1A). S fumigation decreased Chl contents (P < 0.05) in L12 (Figure 1B) and significantly reduced the ΦPSII of leaves (P < 0.05). However, S fumigation slightly but insignificantly decreased the Fv/Fm value (Figure 1C). These results indicated that the effects of S fumigation on the photosynthetic system are greater in L12, compared with L0. Meanwhile, S fumigation increased SOD and POD activities in L0, whereas the activities decreased in L12 (P < 0.05) (Figure 1E). Therefore, the antioxidant enzyme activities in cucumber leaves first increased then decreased with S fumigation. Moreover, the MDA contents in L12 significantly increased (P < 0.05), but L0 only showed slight differences (Figure 1D).

Sulfur fumigation had no significant impact on fruit phenotype (Figure 2A) and on fruit development (Figure 2B). The single fruit mass and the peel color of F12 were not significantly different (Figures 2C,D). Moreover, S fumigation significantly reduced the TSS content in F12 (P < 0.05) (Figure 2E). Overall, S fumigation did not affect cucumber development and appearance but reduced the nutritional value of commercial fruits.

A total of 157.27 Gb cleaned transcriptome reads with 94.74–95.28% of bases scoring Q30 was generated from 24 samples, including three biological replicates per sample (L0, L12, F0, and F12) (Supplementary Table 3). Transcriptome analysis indicated that the leaves had more DEGs than fruits. Most of the DEGs in leaves were upregulated, while most of the DEGs in fruits were downregulated (Figure 3A). The GO enrichment analysis of the DEGs identified in L0, L12, F0, and F12 provided additional gene functions (Figure 3C). The GO terms related to abiotic stress responses in plants, including response to stimulus, response to oxygen-containing compounds, response to hydrogen peroxide, and response to reactive oxygen species (ROS), were enriched. In addition, the secondary metabolite biosynthetic processes, including organic acids and lipids, and hormone signal transduction, especially abscisic acid (ABA), were enriched in the leaves and fruits. These results indicate that S fumigation influenced more transcript expression in leaves than fruits, and the effects were significantly related to abiotic stress responses.

Notably, 34 DEGs were common between leaves and fruits (Figure 3B). Interestingly, expression patterns of 33 (of the 34-common leaf-fruit) DEGs were consistent at F0 and F12. Meanwhile, several plant abiotic stress-related genes were differentially expressed in all the groups. Several genes, including calcium uptake protein 1 homolog, mitochondrial-like (CsaV3_2G003630), octicosapeptide/Phox/Bem1p domain-containing protein kinase (CsaV3_2G024940), serine-rich protein-like protein (CsaV3_2G003620), auxin-repressed protein (CsaV3_6G050280), and NAC domain (CsaV3_3G041280), were downregulated in L0 and upregulated in L12 (Supplementary Table 4). These genes were differentially expressed in L0 and L12, reflecting different stress responses under short-term and long-term S fumigation. More importantly, inositol oxygenase (CsaV3_1G040000) and F-box protein-like (CsaV3_1G038410) were upregulated in all four groups, indicating alleviated oxidative damage (Duan et al., 2012), and induced protein degradation mediated by ubiquitination (Craig and Tyers, 1999).

Notably, 36 WGCNA modules were identified and hierarchically clustered as shown in different colors (Figure 4A). The Pearson’s correlation between modules and phenotypic traits was illustrated in the heat map (Figure 4B). Within the 36 co-expressed gene networks, the green-yellow module (1,630 genes) showed the highest positive correlation with SL0 (r² = 0.78 and P = 7.00E-06). The cyan module (1,835 genes) showed genes significantly associated with SL12 (r² = 0.76 and P = 1.00E-05) (Figure 4B).

The networks of selected key WGCNA-identified genes related to stress response were constructed using Cytoscape (Figures 4C,D). The co-expression network genes with respective annotations are listed in Supplementary Tables 5, 6. Multiple transcription factors, including WRKYs, NACs, ERFs, HSFs, MYBs, and protein kinases (i.e., LECRK3, WAK2, CRK25, CRRSP38, PERK1, CPK32, and XA21) were identified in the green-yellow and cyan modules of SL0 and SL12. Interestingly, transcripts related to phenylpropanoid biosynthesis and calcium signal transduction, including PAL, CCR, and CP, CML, were identified in the green-yellow module. HPL (allene oxide synthase) and LOX (lipoxygenase) related to jasmonic acid biosynthesis were selected in the cyan module. Therefore, S fumigation triggers various defense reactions in cucumber leaves, including phenylpropanoid biosynthesis and jasmonic acid biosynthesis.

The LECRK3, RLK5, poxN1, ASPG1, AMC1, PERK1, XA21, HPL, CYSK, and PHOS32 were identified as hub genes, which included multiple protein kinases and various oxidative stress-responsive genes (Table 1). Gene ASPG1 participates in ABA-mediated stomatal closure and improves stress resistance (Yao et al., 2012), while AMC1 is regulated by calcium signal transduction (Zhu et al., 2020). Meanwhile, poxN1 and XA21 genes mediate various plant defense responses (Kim et al., 2011; Luo et al., 2012). The CYSK gene, related to the synthesis of S-containing amino acid (cysteine), was also identified as a hub gene. Therefore, S fumigation possibly impacts the cysteine metabolism of cucumber. The RT-qPCR expression for the hub genes was consistent with the transcriptome, confirming the reliability of the transcriptome data (Supplementary Figure 1). Overall, various defense regulatory responses such as ROS scavenging, ABA signal transduction, and calcium ion oscillation constituted the cucumber immune defense response after S fumigation.

Previous analyses illustrated that S fumigation significantly affected the ABA signal transduction pathway. Inhibited by ABA receptor PYLs, the protein phosphatase 2Cs (PP2Cs) was downregulated, activating protein kinase SnRKs (Figure 5). Interestingly, the expression of these core factors was consistent in L0 and L12 under S fumigation, indicating that both short-term and long-term S fumigation significantly activated the primary ABA signal transduction pathway. Downstream of the ABA signal transduction pathway and activation of SnRKs upregulated the NADPH oxidase RBOHF profile, which promoted the accumulation of superoxide free radicals. E3-ubiquitin-ligase HOS attenuated the SnRK profile, whereas PUB inhibited the expression of PP2Cs. Furthermore, SnRKs caused differential expression of WRKYs, NACs, and MYBs. These results suggest that protein-level ABA signaling was adapted after S fumigation.

Meanwhile, other downstream defense transcripts of the ABA signal transduction pathway, except HOS, were upregulated in L0 and downregulated in L12 by S fumigation. This means that short-term and long-term S-induced stress responses are different. Moreover, the ABA signal transduction pathway was significantly downregulated in the S fumigated F0 but not F12, possibly indicating that S fumigation affected endogenous ABA metabolism in the ovary.

The metabolic profiles were identified by UPLC-MS/MS. The PCA score plots showed that the L0, L12, and F0 groups were clearly separated from the control along with the first principal component. However, there was no noticeable separation between the F12 group and the control group (Supplementary Figure 2A). To maximize the separation between groups, OPLS-DA was performed on the basis of the metabolite data. The score plots showed that L0, L12, F0, and F12 groups were clearly separated from the control (Supplementary Figure 2B). This indicated that S fumigation altered metabolite profiles in cucumber. The Pearson’s correlation coefficient between each sample was calculated to measure the repeatability and quality (Supplementary Figure 3). The one-way ANOVA results showed that L12 had more differentially expressed metabolites (DEMs) than L0, whereas F12 had fewer DEMs than F0 (Supplementary Figure 4). Overall, there were more downregulated than upregulated DEMs (Supplementary Figure 4A).

The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that the glycerophospholipid metabolism, arachidonic acid metabolism, and linoleic acid metabolism were substantially enriched in leaves and fruits (Figure 6). Additionally, sucrose and starch, galactose, lysine, arginine metabolic pathways, and the TCA cycle were enriched in leaves. Supplementary Tables 6, 7 summarize the significantly changed metabolites mostly related to S fumigation in leaves and fruits. Various sugars, glycosides, organic acids, and amino acids were significantly changed in leaves and fruits (P < 0.05) (Supplementary Tables 7, 8). These results indicate that S fumigation regulates carbon and nitrogen metabolism in cucumbers.

The Pearson’s correlation analysis of the top 50 VIP DEMs in each group identified 21, 9, and 15 lipid metabolites in L0, L12, and F0, respectively. The lipid metabolites positively correlated with sugars (i.e., sucrose, trehalose, galactose, and mannose) (Supplementary Figures 5–7). Interestingly, the lipid metabolites from L0 and F0 positively correlated with several glycosides and negatively correlated with citric acid, maleic acid, and 17-hydroxylinolenic acid (Supplementary Figures 5, 7). However, the lipid metabolites from L12, except PA(15:0/20:3(5Z, 8Z, 11Z)) and TG(16:0/18:1(9Z)/18:1(9Z)), positively correlated with several organic acids [17-hydroxylinolenic acid, pyrrolidone carboxylic acid, and 3,4,5-trihydroxy-6-(2-hydroxy-6-methoxyphenoxy)oxane-2-carboxylic acid] and negatively correlated with several glycosides (Supplementary Figure 6). Thus, short-term and long-term S fumigation had differentially regulated cucumber metabolites.

Citric acid, which participated in the cucumber TCA cycle, was significantly upregulated in cucumber (Figure 7). Furthermore, several sugars, including sucrose, fructose, glucose, raffinose, and galactose, significantly decreased. S fumigation significantly reduced the contents of lipids such as phosphatidic acid, phosphatidylserine, phosphatidylcholine, sphingomyelin, and glycerophosphocholine. More importantly, the contents of triglyceride and proline significantly increased under S fumigation. Triglycerides and proline are essential for peroxide loading (Feussner et al., 1997) and relieving osmotic stress (Abogadallah et al., 2003), respectively. These results proved that S fumigation significantly reduced carbon and lipid metabolism in cucumber, a factor related to the decreased leaf photosynthetic efficiency caused by S fumigation.

To explore the effects of S fumigation on the S-related metabolism of cucumber leaves and fruits, the integration of transcriptome and metabolome was performed using O2PLS. Of note, 57 transcripts and 14 metabolites were used for model construction (R²X = 0.694 and R²Y = 0.838). As shown in Figure 8, CYSK (CsaV3_5G030130), one of the WGCNA-identified hub genes, was positively correlated with S-containing amino acids (i.e., cystine and homocysteine). Meanwhile, strong negative correlations could be seen between glutathione and transcripts involved in sulfate transport, activation, and reduction (i.e., CsaV3_3G045740, CsaV3_1G010690, and CsaV3_6G048960). These genes were significantly downregulated in L12, which indicates that sulfate absorption and assimilation decreased under long-term S fumigation (Figure 9). Moreover, S fumigation significantly increased the contents of various S-containing metabolites, for example, cystathionine, O-phosphohomoserine (OPHS), and S-adenosylmethionine (SAM) in cucumber leaves and fruits. The contents of sulfate, serine, cystine, glutathione, and homocysteine slightly but insignificantly increased. Altogether, these results indicate that the increase of secondary S-containing metabolites in S-fumigated cucumbers may be through the conversion of non-S-containing metabolites, rather than through the enhancement of cucumber sulfate assimilation.