Tomato (S. lycopersicum) ‘Jinpeng No.1’(Ding et al., 2019) was used in this research and experiments were conducted in Chengdu, Sichuan Province, China (104.06°E, 30.67°N) between December 2021 and March 2022. Seeds were sown in a plastic seedling tray (53 × 27.5 × 4.5 cm) filled with substrate (Pindstrup, Demark) and housed within an artificial climate chamber (RTOP-1000D, Top Yun Co. Ltd., Hanzhou, China) with climate settings held at 28 ± 1°C during the day and 18 ± 1°C during the night with 65 ± 5% relative humidity and a photoperiod of 14 h. Three replicated groups, each containing 100 seeds were established. Three weeks after sowing, 60 uniform seedlings from each group with two fully expanded leaves were transplanted into 7 × 7 × 8 cm black plastic pots filled with substrate (Pindstrup, Demark) in a Venlo-type arrangement, with double spans in a north–south orientation greenhouse (9.6×4×4.5m) at a set climate of 30°C/15°C temperatures (day/night) and 65% relative humidity with automatic air conditioning. The daily maximum natural indoor light intensity (PPFD) varied from 100 to 250 µmol·m^-2·s^-1 (measured at the same height of top canopy of tomato plants). Plants were set 10-13 cm apart from one another. After three weeks of irrigation with a half dose of Yamazaki nutrient solution (EC 1.0 ± 0.2 mS/cm), the dose of the solution was doubled (pH 6.5 ± 0.5, EC 2.0 ± 0.5 mS/cm) until the conclusion of the experiment.

120 tomato seedlings grown in the same environment were taken 2 weeks after transplanting and randomly divided into two groups. Light-emitting diodes (LED; Philips Netherlands Ltd.) were used as light sources ( Figure 1A ). The lighting was processed in two orientations: abaxial leafy supplemental lighting (AB) and adaxial leafy supplemental lighting (AD) ( Figure 1B ), according to previous reports (Song et al., 2016; Jiang et al., 2017) with minor adjustments. The LED was fixed on a movable beam to ensure that illumination distance from the adaxial epidermis of the inner canopy truss or the abaxial epidermis of the lowest leaf truss was maintained at 10 cm. For this, the third leaf from both inner canopy truss or lowest truss was taken as a reference, and plant position was adjusted to ensure vertical growth and a consistent plane of stem axis within the same row when necessary. The supplemental lighting PPFD, measured at a distance of 10 cm from the LED module, was 200 µmol·m^-2·s^-1 with a supplemental lighting photocycle of 16h (6:00-22:00), maintained by an integrated digital timer-dimmer-transformer (EEIO-600W-1000W, Shengyuan Electric Appliance Co., Ltd, Zhongshan, China). Each treatment consisted of three rows of plant benches with each row containing 20 plants. After 16 hours of light treatment, the whole leaves of 5 plants from either AB or AD treatment which were randomly selected were mixed as a sample (Liu et al., 2020) for RNA-seq and proteomics analysis.

Gas-exchange measurements were conducted based on previous reports (Song et al., 2016). In brief, we selected the second terminal leaflets of leaves on the fifth youngest node with a portable photosynthesis system (Li-6400XT; Li-Cor Inc., Lincoln, NE, USA) during 11:00–16:00, GMT +8 (9:00–14:00, local time) on the 28th day after transplanting. The net photosynthetic rate (P _N), stomatal conductance (Gs), and transpiration rate (Tr) were measured. Measurements were conducted with PPFD, leaf temperature, CO₂ concentration, and relative humidity of 800 ± 5 µmol·m⁻²·s⁻¹, 28 ± 1°C, 400 ± 2 µmol·m⁻², and 63 ± 2%, respectively.

A total of 6 S. lycopersicum samples, including 3 AB samples and 3 AD samples, were analyzed by RNA-Seq. Total RNAs were extracted from frozen fresh tomato leaves using an EASYspin Plus Kit according to the manufacturer’s instructions (Aidlab Biotechnologies Co. Ltd., Beijing, China). The quality and quantity of extracted RNAs were measured using agar gel electrophoresis and Nanodrop micro spectrophotometry in combination (Thermo Scientific, Wilmington, DE, USA). RNAs from three biological replicates (0.5 g per sample) across at least five plants with the same concentration and volume were combined for RNA-seq. The sequencing library was constructed using a NEBNext Ultra RNA library prep kit (NEB#E7530, New England Biolabs, Ipswich, MA, USA). The quality of the cDNA library was measured using a DNA 1000 assay Kit (5067-1504, Agilent Technologies, Santa Clara, CA, USA) prior to sequencing on an Illumina HiSeq TM 2500 by Gene De novo Biotechnology Co. (Guangzhou, China). RNA-seq data was downloaded from the SRA database (accession number: PRJNA895868). Clean reads were compared to the reference genome sequence using HISAT software (Kim et al., 2015).

Differentially expressed genes (DEGs) between AB samples and AD samples were identified using the DEGseq software package (http://www.bioconductor.org/packages/2.6/bioc/html/DEGseq.html). Manually identified DEGs (log2 value≥1.5-fold difference, p-value less than 0.01) were then subjected to enrichment analysis using Gene Ontology (GO) functions and KEGG pathways. GO DEGs enrichment analysis provided all GO terms that were significantly enriched in DEGs compared to the genomic background. All DEGs were mapped to GO terms in the Gene Ontology database (http://www.geneontology.org/). Significantly enriched GO terms (FDR correction p-value ≤ 0.05) were identified by a hypergeometric test by comparing them to the genomic background. Pathway enrichment analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Pathways with FDR-corrected p-values ≤ 0.05 were defined as significantly enriched DEG pathways.

At least 5 seedlings were mixed in each replicate. Total protein was extracted using the cold acetone method and labeled with iTRAQ tags. Shotgun proteomic analyses were performed using an EASYnLC™1200 UHPLC system (Thermo Fisher, Shanghai, China) with an Orbitrap Q Exactive HF-X mass spectrometer (Thermo Fisher, Shanghai, China). iTRAQ quantification was implemented using IQuant software (Wen et al., 2014). Proteins with a fold change of > 1.2 or < 0.8 and unadjusted significance level p < 0.05 were considered differentially expressed proteins.

The enrichment analysis of differentially expressed proteins (DEPs) was performed by GO and KEGG analysis. The iTRAQ proteomic data were deposited in the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the iProX partner repository, using the iProX data license number is PXD038211.

A set of DEGs and DEPs identified in this research (30 total genes, 25 genes from RNA-seq and 5 from iTRAQ) were selected were selected for qPCR analysis and verification of transcriptional changes after AB treatments for 0, 1, 2, 4, 6, 8, 12, and 24 hours. AD plants were maintained as controls. Tomato actin was used as an internal reference. The primers used were designed using Primer Premier 5.0 (Premier) and are listed in Table S1 .

Triplicate data were analyzed using SAS 9.0 software (SAS Institute Inc., North Carolina, USA) according to SAS Tutorials: Analyzing Data (https://libguides.library.kent.edu/SAS/AnalyzeData). The statistical significance of the difference was evaluated by a Student’s t test and least square means analysis at the level <0.05.

The morphology of tomato seedlings after 2 weeks of the abaxial leafy supplemental lighting (AB) and the adaxial leafy supplemental lighting (AD) can be seen in Figure 2A . Compared with AD, AB significantly increased tomato plant height, while fresh weight increased slightly but not significantly ( Figures 2B, C ). In addition, AB significantly increased the P _N and Gs ( Figures 2D, E ), while Tr was increased slightly but not significantly ( Figures 2F ).

RNA-Seq analysis yielded an average of 11.23 Gb of data per sample. The average alignment rate of the sample comparison genome was 95.17%, and the average alignment rate of the compared gene set was 88.03%. The number of predicted new genes was 25,970 and the total number of detected expressed genes was 47,249, of which 22,279 were known, and 24,970 were predicted new genes. A total of 33,220 new transcripts were detected, of which 503 belonged to novel alternatively spliced isoforms of known protein-coding genes, and 25,970 belonged to new protein-coding gene transcripts. The remaining 6,747 belonged to long non-coding RNAs.

The screening conditions for DEGs were FDR < 0.05 and |log2FC| > 1. A total of 10,998 genes were differentially expressed, with 5280 down-regulated and 5718 up-regulated genes ( Tables S2 ). MA plot, Volcano plot, and Scatter-plot were used to display the distribution of DEGs, and an expression heat map was made for each group of DEGs, shown in Figure S1 .

GO and pathway enrichment analyses were performed for all significantly DEGs ( Figure 3 ). Different comparisons exhibited similar distribution patterns with regard to the numbers and types of enriched pathways, which may be divided into three main functional groups, including 25 biological processes, 16 molecular functions, and 12 cellular components ( Figure 3A ). Significant differences were observed at the level of gene expression in multiple pathways, chiefly including cell process (2564 genes), metabolism process (2480 genes), biological regulation (1017 genes), regulation of biological process (905 genes), response to stimulus (793 genes), and cellular component organization or biogenesis (487 genes).

According to the DEG results, we carried out KEGG biological pathway classification and enrichment analysis. The pathway classification findings illustrated that the functions of DEGs were mainly concentrated in five branches, including Cellular Processes, Environmental Information Processing, Genetic Information Processing, Metabolism, and Organic Systems ( Figure 2B ). The majority of gene functions were focused in metabolic pathways. Pathway enrichment results showed that the top five DEGs were predominantly concentrated in plant hormone signal transduction, cutin, suberine and wax biosynthesis, mismatch repair, sphingolipid metabolism, and other types of O-glycan biosynthesis ( Figure 3C ). Further analysis demonstrated that the DEGs could be classified into 30 categories, and the top five groups were plant hormone signal transduction, MAPK signaling pathway, starch and sucrose metabolism, and phenylpropanoid biosynthesis ( Figure 3D ).

We made predictions for DEGs with the ability to encode transcription factors (TFs), and, classified and counted transcription factor families to which the differently expressed genes belonged. Our findings showed that the six most abundant transcription factors were MYB, MYB-related, bHLH, AP2, MADS, and NAC which contained 486, 380, 274, 264, 169, and 155 transcription factors, respectively ( Figure S2 ).

We also found that many plant hormone signal-regulated genes were significantly upregulated upon AB treatment, including three gibberellin biosynthesis genes 20ox-3, GA2ox5, and GAI, four auxin transport genes PIN4, PIN6, PIN7, and PIN9, two auxin response genes IAA13 and IAA23, five auxin response factors ARF5, ARF8, ARF9, ARF12, and ARF18. Special regulatory genes, such as three sugar transporters SWEET1, SWEET12, and SWEET14, four cell cycle regulators cdkB2, CycA1, CycA2, and cycd3c3, were also significantly upregulated under AB treatment ( Table S2 ).

To investigate the effect of AB on protein expression, we performed iTRAQ analysis on tomato leaves under AB and AD treatments. Across all samples, a total of 338,790 secondary spectra were generated and downloaded. Using the filter standard of “1% FDR”, a total of 21,668 peptides and 5,390 proteins were identified. The significantly different proteins (DEPs) were identified as having a fold change> 1.2 or <0.8 and a Q-value < 0.05, and 152 DEPs, including 65 upregulated proteins and 87 downregulated proteins, were identified ( Table S3 ).

We further conducted functional classification and enrichment analyses of these identified DEPs. The biological processes of DEPs often included cellular process, metabolic process, biological regulation, response to stimulus, regulation of biological process, and cellular component organization or biogenesis ( Figure 4A ), which was consistent with GO analysis results of DEGs ( Table S2 ). Pathway enrichment analysis of DEPs demonstrated that the altered biological pathways were mainly distributed in oxidative phosphorylation, AGE−RAGE signaling pathway in diabetic complications, arginine and proline metabolism, mRNA surveillance pathway, phagosome, and plant hormone signal transduction, et al. ( Figure 4B ).

Integrating proteomic and comprehensive transcriptomic data analysis provided an important validation tool for the expression of key genes. However, it has been established that changes in gene expression do not imply a corresponding change in protein content. We integrated transcriptomic and proteomic data to analyze the regulation of gene expression changes in response to protein expression changes. Our findings indicated that 17 genes and proteins were upregulated together, including NADH dehydrogenase, aspartate aminotransferase 3 (ASP3), casein kinase II subunit alpha, stromal cell-derived factor 2-like protein (SDF), protein argonaute 5 (AGO5), xyloglucan endotransglucosylase/hydrolase 1(XTH1), putative ABC1 protein, metal transporter Nramp3 (NRAMP3), PLAT domain-containing protein 1 (PLAT1), lysophospholipase BODYGUARD 3 (BDG3), formate dehydrogenase (FDH1), thiamine pyrophosphokinase 1 (TPK1), inositol phosphorylceramide glucuronosyl transferase 1 (IPUT1), patatin-like protein 1(PLP1), LysM domain receptor-like kinase 3 (LYK3), Non-specific lipid-transfer protein 2 (LE16), and heat stress transcription factor A-1 (HSFA1) ( Table S3 ).

Moreover, 18 genes and proteins were downregulated simultaneously, including phenylalanine ammonia-lyase (PAL5), pre-mRNA-splicing factor ATP-dependent RNA helicase (DEAH4), histone H2A.1, serine/threonine-protein kinase (SAPK3), DNA-directed RNA polymerases II, IV, and V subunit 9A (NRPB9A), KH domain-containing protein, protein SMAX1-LIKE 3 (SMXL3), 40S ribosomal protein S16 (RPS16), receptor-like protein kinase (HSL1), polyprenol reductase 2 (PPRD2), basic blue protein, nuclear pore complex protein (NUP1), eukaryotic translation initiation factor 3 subunit H (TIF3H1), histone H2A.1, Chlorophyllase-2 (CLH2), exocyst complex component (EXO70A1), and transcription termination factor (MTERF8).

We used qPCR to further determine gene expression levels in S. lycopersicum under AB treatment. The results were consistent with those from RNA-seq and proteomic analyses ( Figure 5 ). It was shown that AB does promote upregulated expression of some key genes, such as 20ox-3, GA2ox5, GA2ox8, GA2ox10, GAI, IAA13, IAA23, PIN4, PIN6, PIN7, PIN9, SWEET12, SWEET14, CycA1, CycA2, and Cryptochrome DASH. Our results ultimately contribute to the understanding of plant light response mechanisms and the discovery of efficient light utilization genes.