Comparative physiological and transcriptomic analysis of pear leaves under distinct training systems

Canopy architecture is critical in determining the light interception and distribution, and subsequently the photosynthetic efficiency and productivity. However, the physiological responses and molecular mechanisms by which pear canopy architectural traits impact on photosynthesis remain poorly understood. Here, physiological investigations coupled with comparative transcriptomic analyses were performed in pear leaves under distinct training systems. Compared with traditional freestanding system, flat-type trellis system (DP) showed higher net photosynthetic rate (PN) levels at the most time points throughout the entire monitored period, especially for the interior of the canopy in sunny side. Gene ontology analysis revealed that photosynthesis, carbohydrate derivative catabolic process and fatty acid metabolic process were over-represented in leaves of DP system with open-canopy characteristics. Weighted gene co-expression network analysis uncovered a significant network module positive correlated with PN value. The hub genes (PpFKF1 and PpPRR5) of the module were enriched in circadian rhythm pathway, suggesting a functional role for circadian clock genes in mediating photosynthetic performance under distinct training systems. These results draw a link between pear photosynthetic response and specific canopy architectural traits, and highlight light harvesting and circadian clock network as potential targets for the input signals from the fluctuating light availability under distinct training systems.


Materials and methods
Plant materials. Ten-year-old 'Wonhwang' (Pyrus pyrifolia Nakai cv. Wonhwang) pear trees growing in the experimental orchard (30.292°N, 114.143°E) of Hubei Academy of Agricultural Sciences, Research Institute of Fruit and Tea, were used for this study. The planting distance was 3 m between the rows and 4 m between the trees. The trees in this experimental orchard received fertilizers, irrigation, and chemical thinners in accordance with the local recommendations. Trees have been trained to SP or DP system as described previously 9 . Briefly, the SP system consists of one central stem and 3-5 primary branches that are upright and located at about 0.5-2 m from the ground. Both the central stem and primary branches have several smaller sub-branches. The DP trellis system is near to a Y-shaped system, which features a support structure consisting of one central stem and two primary branches bent in opposite directions along the row. The height of the trellis is 1.7-1.8 m. At planting, trees are headed at 1.2-1.3 m above the ground so that they would have an upright central leader at the base. Heading back produced many strong shoots that were selected and trained to the two proper arms of the trellis system. Each primary branch with an incline of 45° above the horizontal has equally spaced subbranches that were naturally tied on horizontal steel wires suspended from concrete posts. The experiment was a randomized complete block design with three replications. Trees within each block were randomly selected, which represented biological replicates per training system. Each tested tree was divided into eight leaf locations, i.e., sunny side ( Photosynthetic measurements. Diurnal courses of photosynthetic parameters including P N and leaf temperature (T) were measured using the portable TPS-2 photosynthesis system (PP. Systems Inc., USA). All measurements were carried out at regular intervals of 2 h (n = 5), between 08:00 and 16:00 on sunny and clear days, during the spring-summer productive period including 15 DAF (day after flowering), 45 DAF, 75 DAF and 105 DAF. Leaves on the East side at 8:00, 10:00, 12:00 and leaves on the West side at 14:00, 16:00 are considered Scientific Reports | (2020) 10:18892 | https://doi.org/10.1038/s41598-020-75794-z www.nature.com/scientificreports/ as sun leaves, while leaves on the West side at 8:00, 10:00, 12:00 and leaves on the East side at 14:00, 16:00 are considered as shady leaves. Photosynthetic photon flux density (PPFD) was measured with LI-180 Spectrometer (LI-COR Inc., USA) at 5 cm above the surface of the leaves in specific location. All PPFD measurements were taken every 2 h between 08:00 and 16:00 on sunny and clear day (105 DAF). For each biological replicate of each leaf location, photosynthetic parameters were measured on three leaves (three technical replicates), which were youngest, fully expanded sixth leaves from 1-year-old shoots, at five time points of the day. Statistical analyses were performed with a model of one-way ANOVA analyses (IBM SPSS Statistics 19 software), followed by the Duncan's multiple range tests at p < 0.01 as described previously 32 .
RNA isolation and RNA-seq. Leaf 34 . Correlation analysis for gene expression levels between any two samples was performed by spearman algorithm. Differentially expressed genes (DEGs) were detected using DESeq2 35 , and Benjamini-Hochberg's method was used to control the false discovery rate (FDR). Only genes with absolute value of log 2 (fold change) ≥ 1 and adjusted p value < 0.05 were considered significantly differentially expressed. We analyzed the gene relationships and identified the overlapping DEGs using VennDiagram 36 .

Quantitative real time-PCR (qRT-PCR) analysis.
To validate the RNA-seq results, expression of selected genes was determined by qRT-PCR analysis. One microgram of total RNA was reverse-transcribed using RevertAid First Strand cDNA Synthesis Kit (Fermentas, USA). qRT-PCR primer pairs specific to selected genes were designed by Primer Premier 5.0 software (Supplementary Table S1). The primers were further confirmed with a melting curve analysis after amplification of each tested genes. Two reference genes, i.e. PpSKD1 and PpYLS8, which proved to be stably expressed in leaves during different growth stages under distinct pear training systems, were used as suitable internal controls to normalize the qRT-PCR data 9 . Relative quantification was calculated based on the Ct method (2 −△△CT ). Three independent biological replicates for each sample and three technical replicates of each biological replicate were performed.

Functional annotation, pathways and transcription factors (TFs) analysis.
In order to identify the significantly enriched Gene Ontology (GO) terms, each of the DEG lists were mapped into GO terms in the database (https ://www.geneo ntolo gy.org/). Go functional enrichment analysis was performed using Goatools software by Fisher's exact test 37 . Terms are considered enriched at p < 0.001. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched analysis were performed (https ://www.genom e.jp/kegg/), using the criterion of a Benjamini-Hochberg corrected p value < 0.05 [38][39][40] . Up-regulated and down-regulated TFs were identified and classified into different families using PlantTFDB v4.0 (https ://plant tfdb.cbi.pku.edu.cn), with a threshold E-value < 1E−5.
Weighted gene co-expression network analysis (WGCNA). A network analysis based on gene-totrait correlations was performed using WGCNA R package 41 . As it is believed that low expressed and non-changing genes provide limited information in a co-expression network building, 13,521 pear genes were selected based on their expression values ≥ 1 FPKM in one or more samples and coefficient of variation (CV) ≥ 0.1. A soft threshold value, power of 6, was used to transform the adjacency matrix to meet the scale-free topology criteria for optimal clustering. Modules whose eigengenes were highly correlated were merged with a mergeCutHeight of 0.25. The minimal module size was set to 30 genes. Module-trait associations were assessed by calculating the spearman's correlations between the module eigengenes (MEs) and the P N values. Top 30 genes in intramodular connectivity ranking were counted as hub genes in a given module. The network for hub genes was visualized using the Highcharts (v6.0.4) (https ://www.highc harts .com/). After completing the P N measurements, the same leaves were excised to measure leaf mass per area. The leaves were then oven-dried at 60 °C for one week to determine their dry mass, and the samples were digested with  42 . Soluble protein content of leave was measured by the Coomassie brilliant blue G250 method. Leaf sample was homogenized in 1 ml distilled water, and the solution was centrifuged at the speed of 10,000 rpm for 10 min at 4 °C. The supernatant was used for colorimetric assay of soluble protein content at a wavelength of 595 nm. Standard curve were prepared with bovine serum albumin (Sigma, ultra 99%).

Measurement of photosynthetic pigments
PpRCA enzyme activity was determined by an enzyme-linked immunosorbent assay (ELISA) using Rubisco Activase Assay Kit according to the manufacturer instructions. The absorbance was measured at 450 nm. The activity of RCA enzyme in samples was calculated by the standard curve.

Results
Spatial and temporal variability of biochemical photosynthetic parameters under distinct training systems. To understand the response of photosynthesis to variable canopy positions, the P N fluctuations of eight typical canopy locations were monitored from 15 to 105 DAF. The leaf P N varied between 1.30 and 20.27 µmol m 2 s −1 during the time period (Fig. 1A). Among these, the sun leaves of DP at 45 DAF (DPSU45) showed the highest daily median value (16.57 µmol m 2 s −1 ), whereas the shade leaves of SP at 105 DAF (SPSH105) presented the lowest average P N (1.97 µmol m 2 s −1 ). The difference may due to the sensitivity to light depending on the leaf growth stages and spatio-temporal variations in light exposure.
The leaf P N in sunny side displayed relatively higher variation, with P N values 2.00-20.27 µmol m 2 s −1 (Fig. 1A). In contrast, the leaf P N in shady side showed relatively narrow range spanning 1.30-13.47 µmol m 2 s −1 . The rates of P N were markedly higher in sun leaves than in shade leaves ( Fig. 1B), which is probably due to higher light intensity in sunny side.
Fruiting leaves and vegetative counterparts had similar photosynthetic capacity at the most time points throughout the entire monitored period, with a few minor exceptions that P N value of SU-EX exhibited higher levels in vegetative leaves than that in fruiting leaves at 75 DAF (SP) and 105 DAF (Fig. 1B).
To analysis the diurnal profiles of photosynthetic activity, the leaf temperature, PPFD and P N parameters were monitored along with five time points throughout daylight ( Fig. 2 Figure S1). Similar changing patterns and levels of leaf temperature were observed between sunny side and shady side (Supplementary Figure S1A). Leaf temperature was increased rapidly in the morning and showed a relatively high peak value at 12:00 h or 14:00 h, and then slight decreased thereafter. The diurnal patterns of the PPFD also showed a variation similar to the changes of leaf temperature, which exhibited its daily maximum at solar noon, except for SP leaves in the interior part of the canopy, where PPFD values remained relatively constant during the course of the day (Supplementary Figure S1B). When looking at the discrimination between the distinct training systems, significantly higher PPFD values were observed in DP leaves than in SP leaves under SU/SH-IN and SH-EX conditions. Under the fluctuating light intensity and temperature, the diurnal patterns of P N showed a strong midday depression in most sun leaves ( Fig. 2A). In shady side, the diurnal course of P N showed approximately 'dome shaped' pattern ( Fig. 2B). Most of them reached a maximum at noon, whereas others displayed no significant variation trend.

and Supplementary
For DP system, leaf P N was not statistically different between exterior and interior part of trees at the most canopy positions (Fig. 1B). In contrast, in sunny side, the SP system displayed significantly higher values in exterior part of the canopy except at 15 DAF stage. More importantly, P N showed relatively higher levels in DP compared with SP, especially for the interior part of the canopy (Figs. 1A and 2), indicating that mutual shading may be more severe in SP system. The difference in P N between SP and DP training systems appeared at the 12:00 h of 15 DAF but became more evident at the later stages (45 DAF, 75 DAF and 105 DAF) ( Fig. 2A). It is conceivable that the DEGs between SP and DP leaves in SU-IN at the last three stages may play important roles in determining the differences of photosynthetic efficiency under distinct training systems; hence, these were subject to detailed investigation.
Overview of the RNA-Seq data. To explore how light availability affected by distinct canopy structures impacts transcription at the genomic scale, we performed RNA-seq analysis comparing gene expression changes between DP and SP leaves in SU-IN at three successive stages (45 DAF, 75 DAF and 105 DAF) of leaf development. After removing low quality reads, a total of 41.85-53.26 million clean reads were generated from each library (Supplementary Table S2). The percentages of Q20 were higher than 98.3% in all libraries, representing high quality sequencing. We found 74.49-78.38% of the clean reads in the libraries were mapped onto the pear reference genome. Based on FPKM value, we determined the number of genes expressed (FPKM > 1) in individual leaf sample (Fig. 3A). In total, 24,734, 24,203 and 24,886 genes were found to be expressed in SP leaves at 45 DAF, 75 DAF and 105 DAF, respectively. Similarly, 25,524, 24,891 and 24,937 genes were identified in the samples from the respective stages of DP leaves. We also observed similar distributions of gene expression levels across all samples. Approximately 55.7% of expressed genes were in the 1-10 FPKM range, and 39.2% of expressed genes were in the range 10-100 FPKM.
To investigate the relationships among the samples used in this study, we performed hierarchical clustering based on the whole-gene expression data set. The result showed that all samples had a high degree of similarity (more than 0.883) and clustered in three discrete groups (Fig. 3B). This analysis clearly revealed the distinctness of the early stage ( www.nature.com/scientificreports/ presence of different transcriptional programs. In addition, among the two later stages samples, the SP system leaf samples (SP75A/B/C and SP105A/B/C) clustered closely. Likewise, DP system leaf samples (DP75A/B/C and DP105A/B/C) also showed very tight clustering. This indicates distinct transcriptional programs active between the SP and DP samples, but somewhat similar within each other. All the stages of leaf development analyzed in this transcriptome study showed significantly different photosynthetic rates, which imply dramatic differences in their transcriptional programs. Three pairwise transcriptome comparisons were performed between the SP and DP leaves. A total of 203 (168 up-regulated and 35 down-regulated), 3460 (1603 up-regulated and 1857 down-regulated) and 1358 (431 up-regulated and 927 down-regulated) genes were significantly differentially expressed at 45 DAF, 75 DAF and 105 DAF, respectively (Fig. 3C). Consistent with correlation analysis results, striking differential expression was observed at the two later stages, which may be due to distinct canopy architectures present significantly different microclimate during that growth stages.
Venn diagram analysis showed that only 11 genes were commonly up-regulated in the three stages (Fig. 3D). There was no down-regulated gene overlapped among all the three stages, but 4 and 692 genes were commonly down-regulated in 45 DAF/75 DAF pair and 75 DAF/105 DAF pair, respectively. More overlapped genes were observed between 75 and 105 DAF than between 45 and 75 DAF in both up-regulated and down-regulated genes. Asterisks below indicate significantly higher levels of P N between SP and DP system. www.nature.com/scientificreports/ The expression levels of nine interesting DEGs (five up-regulated and four down-regulated) were validated by qRT-PCR (Fig. 3E), all of which are known to be related to electron transport (PpPPO, polyphenol oxidase, LOC103927471), carbohydrate metabolic process (PpG3PDH, glycerol-3-phosphate dehydrogenase, LOC103927132 and LOC103953210), photorespiration (PpRbcS, the small subunit of ribulose-1,5-bisphosphate carboxylase, LOC103948502), as well as some TFs (PpPIF4, putative bHLH transcription factor, LOC103947396; PpDOF, putative DNA-binding one zinc finger transcription factor, LOC103927195; PpMYB4, putative MYB transcription factor, LOC103931820; PpERF, putative ethylene responsive factor, LOC103952763; PpHSF, putative heat stress transcription factor, LOC103960239). qRT-PCR expression profiles of the selected genes were consistent with the deep sequencing results, indicating that the RNA-seq data in this study is reliable.
Analysis of differentially expressed genes (DEGs) responding to microclimate conditions under distinct training systems.. Go enrichment analysis revealed that ten biological processes were enriched in up-regulated genes in DP leaves at 45 DAF (Fig. 4A and Supplementary Table S3). The up-regulated genes were mainly associated with 'phenylpropanoid metabolic process' (GO:0009698), 'oxidation-reduction process' (GO:0055114) and 'carbohydrate derivative catabolic process' (GO:1901136). Special interest child terms of these GO terms were 'lignin catabolic process' (GO:0046274) and 'glycerol-3-phosphate catabolic process' (GO:0046168) (Fig. 4A,D). Two PpG3PDHs (LOC103927132 and LOC103953210) encoding glycerol-3-phosphate dehydrogenase, two PpCesAs (LOC103958831 and LOC103959701) encoding cellulose synthase A and three PpLACs (LOC103959581, LOC103949457 and LOC103930925) encoding laccase all contribute to enrichment of these GO terms. Thus, canopy structure of DP appears to establish a microclimate favorable to carbohydrate metabolism, in comparison to SP counterpart. While six GO terms were enriched for down-regulated genes, and they are mostly associated with 'sulfate reduction' (Go:0019419) and 'cell redox homeostasis' (GO:0045454) (Supplementary Table S3).
Further, we analyzed the enrichment of GO biological terms at 75 DAF stage. The over-represented categories with the highest confidence among the up-regulated genes in the DP leaves were associated with photosynthesis process and nutrient contribution, including 'photosynthesis, light harvesting' (Go:0009765), 'photosynthesis, dark reaction' (GO:0019685), 'fatty acid biosynthetic process' (GO:0006633) and 'protein transport' (GO:0015031) (Fig. 4B,D). It is worth noting that as many as 17 genes encoding light-harvesting chlorophyll a/b-binding (LHC) protein were clustered in the DP leaves enriched category 'photosynthesis, light harvesting' . On the other hand, 65 biological processes were significantly over-represented in the down-regulated DEGs including 'nucleic acid metabolic process' (GO:0090304), 'macromolecule biosynthetic process' (GO:0009059), 'phosphorus metabolic process' (GO:0006793) and 'serine family amino acid metabolic process' (GO:0009069) (Supplementary Table S3). This GO enrichment analysis revealed DP canopy microclimate possibly promotes or triggers the light harvesting response and accelerate fatty acid biosynthesis, while the microclimate conditions of SP canopy architecture activate metabolic pathways.
In addition, six and 16 biological processes were highly enriched among the up-regulated and down-regulated DEGs in the DP leaves at 105 DAF, respectively ( Fig. 4C and Supplementary Table S3). In the group of up-regulated DEGs, one interesting Go category, 'pigment metabolic process' (GO:0042440), was significantly enriched in the DP leaves. Within the down-regulated genes, 'flavonoid biosynthetic process' (GO:0009813), 'nucleic acid metabolic process' (GO:0090304), 'l-ascorbic acid biosynthetic process' (GO:0019853) and 'brassinosteroid homeostasis' (GO:0010268), were significantly enriched categories in the SP leaves.
To identify biological pathways regulated by the distinct training systems, the KEGG enrichment analysis was performed. Our results showed that a total of six and nine pathways were significant enriched in DP leaves at 45 DAF and 75 DAF stages, respectively, but no up-regulated genes were significantly enriched in any pathway at 105 DAF stage (Supplementary Table S4). These enriched pathways included 'amino sugar and nucleotide sugar metabolism' (ko00520), 'photosynthesis' (ko00195), 'photosynthesis-antenna proteins' (ko00196), 'carbon fixation in photosynthetic organisms' (ko00710) and 'fatty acid elongation' (ko00062). In the other hand, the down-regulated genes were grouped into two, two and five significantly enriched pathways at the three stages, respectively. These pathways included 'circadian rhythm-plant' (ko04712) and 'Plant-pathogen interaction' (ko04626) in our study.
WGCNA analysis reveals potential genes associated with photosynthetic performance under distinct training systems. To further investigate candidate genes related to photosynthetic performance under distinct training systems, we performed a WGCNA analysis with all expression genes. A total of 14 distinct modules with sizes ranging from 95 to 2175 genes were identified (Fig. 5A,B). Analysis of the moduletrait relationships revealed that the purple module comprising 214 genes showed strong positive correlation with P N value of fruiting leaves (r = 0.744, p = 0.0004) and vegetative counterparts (r = 0.649, p = 0.00357), indicating that expression of genes in this module may be related to photosynthetic performance. No significantly enriched KEGG pathway in purple module can be observed but six pathways showed near-enriched (p value < 0.05) (Fig. 5C) Effects of canopy conditions on photosynthetic pigments, PNUE, soluble protein and PpRCA enzyme activity. GO analysis indicated that the canopy structure of flat-type trellis system represents an effective strategy for increasing expression of genes that potentially contribute to 'photosynthesis, light harvesting' (Go:0009765), 'photosynthesis, dark reaction' (GO:0019685), 'carbohydrate derivative catabolic process' (GO:1901136), 'protein localization' (GO:0008104) and 'pigment metabolic process' (GO:0042440) (Fig. 4). To examine the effects of canopy conditions on photosynthetic physiology and dry matter accumulation in the field canopies, we investigated photosynthetic pigments, PNUE levels, soluble protein, and PpRCA enzyme activity between distinct training systems (Fig. 6). Compared to SP system, the activity level of PpRCA enzyme was significantly higher in the leaves of DP system. Although no significant differences were observed, the contents of Chl a and carotenoids, PNUE value, soluble protein content, the ratio of Chl a to Chl b (Chl a/b) tended to be higher in DP system, while the content of Chl b showed an opposite trend.

Discussion
Canopy architecture is critical in determining the light distribution, and therefore the photosynthetic productivity of crop. However, information about the mechanism of impact of training systems on photosynthesis is limited. In this study, both GO and KEGG analysis indicated that the up-regulated genes in DP training system participate in photosynthesis processes such as light harvesting (Fig. 4, Supplementary Table S3 and Table S4). The over-representation of genes in these categories was an expected result because of the higher photosynthetic levels in DP systems (Fig. 2). Light harvesting is the first step in the process of photosynthesis, which is regulated by physiological status and environmental signals 43 . LHC proteins principal roles are to efficiently collect light energy and provide photoprotection, which associated with both short-and long-term adaptations to changing environments, such as fluctuating light intensity, temperature changes, or nutrient availability [43][44][45][46] . Some studies indicated that expression levels of LHC genes were induced in response to high light, suggesting a photoprotective function for these genes products 47,48 . Compared with traditional freestanding system with relatively shaded conditions, the flat-type trellis system allows trees to receive sunlight uniformly and opulently 9 .
In the natural environment, DP leaves may suffer excess light that may damage their photosystems, especially in summer (75 DAF-105 DAF). In the present study, 17 PpLhc mRNAs accumulated in the DP leaves with abundant light conditions (Fig. 4D), suggesting a more significant photoprotection and/or light capture roles for DP leaves by enhancing its light-harvesting potential. In higher plants, most photosynthetic pigments (Chls and carotenoids) are coordinated by LHC proteins 49 . Rapid qualitative adjustment in leaf photosynthetic apparatus and photosynthetic capacity can occur in response to short alterations in light availability 50 . As a major change, sun-exposed plants tends to increase the Chl a/b ratio 51 . This general finding was also confirmed in the current study, as a higher photosynthetic levels and Chl a/b ratio were observed in DP leaves with relatively higher irradiance (Figs. 2 and 6B and Supplementary Figure S1B). All these data support that canopy structure of DP is well proportioned and maximize light capture to improve photosynthetic efficiency, which may accelerate nutrient accumulation and transport. Improving the performance of the carbon fixation has the potential to significantly enhance photosynthetic efficiency and crop productivity 52,53 . It is documented that increasing PNUE of the leaves may improve dry matter accumulation 54 . Our results showed that DP system had higher PNUE values and soluble protein contents (Fig. 6C,D), indicating that carbon fixation process may be activated in DP system. Rubisco, which is composed if the eight large subunits (RbcLs) and eight small subunits (RbcSs), is the key enzyme catalyzing the first step of photosynthetic carbon assimilation through the dark reaction 55,56 . Of these, RbcS is thought to be indirectly involved in the catalytic reaction 56 . Overexpression of RbcS genes showed to enhance the catalytic performance of Rubisco, suggesting that RbcS could be an important factor in determining kinetic properties of Rubisco 55,57 . In this study, the GO term 'photosynthesis, dark reaction' was also enriched in DP leaves (Fig. 4B). Three PpRbcS genes (LOC103948502, LOC103948503 and LOC1039346 23) were found to be up-regulated in the DP leaves , which may greatly enhance the catalytic performance of Rubisco and contribute to increase dark reaction activity. We also provide evidence that the activate level of PpRCA enzyme was significantly higher in DP system (Fig. 6E). It is long known that the activation state of Rubisco is regulated by RCA via the maintenance of Rubisco catalytic sites 58,59 . Consequently, PpRCA activity might play an important role in the regulation of PpRbcS genes and photosynthetic carbon assimilation.
In DP system, we noticed a significantly increased expression of genes that were involved in primary metabolism such as carbohydrate derivative catabololic process (GO:1901136), fatty acid and monocarboxylic acid biosynthesis (GO: 0006633 and 0072330) and protein transport (GO:0015031) (Fig. 4), which indicated a more adequate carbohydrate accumulation and utilization for DP system. Though both the leaves in the interior part of the canopy were investigated in the two training system, DP leaves seem to expose to much more irradiance. As compared to shady leaves, sun leaves are often characterized by a higher cutin, lipid and starch content per dry weight basis, and soluble carbohydrate levels 60 . In this study, PpKCS genes encoding 3-ketoacyl-CoA  and LOC103947523) which catalyzes the first reaction in fatty acid elongation and plays a role in cutin and wax biosynthesis 61,62 , were higher expressed in DP leaves (Fig. 4D). Lipids are the major constituents of all membranous structures in plants 63 . The abundant thylakoid membranes in the chloroplast are the site of the photosynthesis light reactions 64 . Lipid biosynthesis genes PpG3PDH (LOC103953210 and LOC103927132) which catalyzes the formation of the backbone of membrane lipids also exhibited higher expression levels in DP system leaves (Figs. 3E and 4D). Overexpression of AtG3PDH gene showed significant enhancement of plastidic lipid content and photosynthetic efficiency 63 . These results imply that enhancement of PpKCS and PpG3PDH activity in DP leaves may promote to lipid accumulation that greatly influences photosynthetic efficiency. Certain stress-assoiated transcription factor families, such as WRKY and NAC genes, have been demonstrated to participate in variable light conditions 65,66 . For example, the AtWRKY22 was found to participate in the darkinduced senescence signal pathway 67 . The ectopic overexpression of TaWRKY7 in Arabidopsis significantly promoted early leaf senescence under darkness treatment 68 . Whole-genome ATH1 Genome Array studies showed that more than one quarter of NAC proteins in Arabidopsis leaves was up-regulated under dark treatment 69 . In accordance to with these observations, increased expression levels of 29 PpWRKY and 22 PpNAC TFs were showed in SP system with lower light intensity compared with the DP system (Supplementary Figure S2). These results implied that senescence-associated TFs, the WRKY and NAC superfamilies, may play crucial roles in regulating, acclimating, and modulating gene expression in photosynthesis process in response to low light. In the other hand, when shaded by neighbouring vegetation, plants are exposed to a variety of light quality signal 70 . PIF TFs have been demonstrated as positive regulators of shade avoidance 71 . The expression levels of two PIF TFs, PpPIF4 (LOC103947396) and PpPIF7 (LOC103931008), were significantly induced in DP system with opencanopy characteristics (Fig. 3E and Supplementary Figure S2). It may suggest that a possible feedback mechanisms mediated by complicated light signaling network may exist due to highly dynamic canopy microclimate.
A new insight into the application of WGCNA and the enrichment of the key module led to the identification of a central role for circadian rhythm regulation related to photosynthetic performance under distinct training system (Fig. 5). Light availability varies not only across the course of the day but also with the changes of canopy architecture, which may provide complex signals to the circadian clock 72 . Recent studies have discovered a link between altered circadian clock regulation and increased levels of photosynthetic activities and biomass in plants [73][74][75] . For example, maize circadian clock genes ZmCCA1s are diurnally up-regulated in the hybrids and target thousands of output genes early in the morning-phased genes; consequently the bindings to carbon fixation genes promote photosynthesis and growth vigor 75 . In our study, network construction highlighted hub genes PpPRR5 and PpFKF1, which showed positive correlation with photosynthetic performance under distinct www.nature.com/scientificreports/ training systems, predicted to play important roles in circadian rhythm (Fig. 5C,D). Arabidopsis homologue AtPRR5, together with other core circadian clock components, form complicated transcription feedback loops that mediating a number of circadian outputs 76,77 . "Red or far-red light signaling pathway" was reported to be an enriched category in AtPRR5 direct-targets 78 . AtFKF1 is a blue light photoreceptor and could target AtPRR5 for degradation 79 . In nature, different canopy conditions combine changes in light spectral quality 50 . The red: far red and blue levels of daylight decreases as it passes through the vegetative canopy 80,81 . In our study, the expression levels of PpPRR5 and PpFKF1 were up-regulated in SP system with mutual shading characteristics. It would be reasonable to speculate that the varied light signal are pronounced at highly dynamic and heterogeneous light conditions, and thus energetically activate the temporal regulation of circadian clock. The oscillator components of the clock, PpPRR5 and PpFKF1, may modulate pear leaves response to light cues and cause a cascade of effects on downstream regulatory pathways, leading to photosynthetic variability.

Data availability
All materials and data sets represented in the current study are available in the main text or the supplementary materials. RNA-Seq data are deposited in the Sequence Read Archive (SRA) database (the accession number PRJNA579772) at the National Centre for Biotechnology Information (NCBI). The metadata are available at https ://www.ncbi.nlm.nih.gov/biopr oject /PRJNA 57977 2.