Transcriptional regulation of photosynthesis under heat stress in poplar

Background Photosynthesis has been recognized as a complicated process that is modulated through the intricate regulating network at transcriptional level. However, its underlying mechanism at molecular level under heat stress remains to be understood. Analysis of the adaptive response and regulatory networks of trees to heat stress will expand our understanding of thermostability in perennial plants. In this study, we used a multi-gene network to investigate the regulatory pathway under heat stress, as constructed by a multifaceted approach of combining time-course RNA-seq, regulatory motif enrichment, and expression-trait association analysis. Results By analyzing changes in the transcriptome under heat stress, we identied 77 key photosynthetic genes, of which 97.4% (75 genes) were down-regulated, and these results conformed to the decreased photosynthesis measured values. According to analysis of regulating motif enrichment, these 77 differentially expressed genes (DEGs) had common vital light-responsive elements involved in photosynthesis. When integrating all the differential expressed genes, 5 co-expressed gene modules (1,548 genes) were identied to be signicantly correlated with 4 photosynthesis-related traits. Thus, based on this, a three-layered gene regulatory network (GRN) was established, which had included 77 photosynthetic genes (in the bottom layer), 40 TFs/miRNAs (in the second layer), as well as 20 TFs/miRNAs (in the top layer), using a backward elimination random forest (BWERF) algorithm. Importantly, 6 miRNAs and 4 TFs were found to be key regulators in this regulatory pathway, emphasizing the signicant roles of TFs/miRNAs in affecting photosynthetic traits. The results imply a functional role for these key genes in mediating photosynthesis under heat stress, demonstrating the potential of combining time-course transcriptome-based regulatory pathway construction, cis-elements enrichment analysis, and expression-trait association approaches to dissect complex genetic networks. Conclusions The heat-responsive pathway in regulating photosynthesis is a multi-layered complex network which is co-controlled by TFs and miRNAs. Our work not only imply a functional role for these key genes in mediating photosynthesis responding to abiotic stress in poplar, but demonstrate time-course transcriptome-based

by the Intergovernmental Panel on Climatic Change report, in the end of 21 st century the climate on the Earth will be warmed by 2-4°C [16]. Therefore, a comprehensive understanding towards relevant gene expression and photosynthetic physiology in response to heat stress contributes to enhancing plant thermostability and limiting the unfavorable impact of climate changes on the productivity of woody crops.
There are numerous studies examining the role of stress in photosystems, electron transport system, photosynthesis-related enzyme activities, pigments, chlorophyll uorescence and gas exchange in plant [17][18][19]. However, a majority of those studies concentrate on plant responses to adapt to the heat stress, while the plant regulatory network in the presence of constant heat stress is rarely mentioned. Populus, a kind of rapidly growing woody plant species, is extensively utilized in wood production and landscaping; but its ecological bene ts could be greatly impacted by increased temperature. Populus have evolved speci c physiological mechanisms for adaption to changes in natural environmental conditions [20,21].
Analyzing the plant regulatory networks and adaptive responses in the presence of heat stress can broaden our knowledge about perennial plant thermostability.
The majority of adaptive response functions by means of gene expression regulation to some extent; as a result, transcript factors responding to heat stress may exert vital parts in the tolerance to abiotic stress [2,22]. Numerous genes can interact mutually as well as with environment, which play certain roles in the heat stress response [3]. A large array of abiotic stress responsive genes has been identi ed in plants with microarray and large-scale transcriptome analyses [23]. These genes play a role in cell protection from damage through producing vital metabolic proteins and enzymes, and in regulating gene expression and signal transduction in the presence of stress [2,3,24]. Of regulating proteins, transcriptional factors (TFs) have critical parts in converting stress signal perception into expression of stress responding genes, which is achieved through the interaction with cis-acting elements different speci c stress responding gene promoter regions during signal transduction. Therefore, it can activate signaling cascade in enhancing plant tolerance to harsh environmental conditions [22,25,26]. About 7% coding sequences in plant genome can be attributed to TFs, a majority of which are in large gene families in comparison with those in yeasts and animals, like the heat stress transcriptional factors (HSFs) family [2,22]. It remains unknown about the numbers and expression of genes involved in trees' photosynthesis regulation under heat stress. As a result, it is of vital importance to characterize and identify genes taking part in plant heat stress tolerance.
Recently, microRNAs (miRNAs), the endogenous short non-coding RNAs, are indicated to exert key parts in the post-transcriptional response to stress in plants [27,28]. MiRNAs function by the negative regulation of gene expression through enhancing target mRNAs degradation or restraining translation. Numerous miRNAs are found to regulate the plant response to both abiotic and biotic stresses, such as nutrient starvation [29], pathogen invasion [30], cold, salt, and drought [31], and oxidative stresses [32]. Although, miRNAs in heat tolerance studies has been made in plants, a genome-wide identi cation of miRNAs responsive to consistent heat stress in poplar has not been conducted. Previous studies on poplar have indicated that many miRNAs are differentially expressed [20,21]. Previous reports on plants under abiotic stress have revealed diverse functions of miRNAs [3,11,33]. For example, through regulation of the target genes, rice TEOSINTE BRANCHED1 and a MYB transcription factor, miR139 is up-regulated and involved in process of rice seedling [28]. Similarly, miR393 was found to participate in responding to drought stress in rice via negative regulation of mRNAs encoding the F-box auxin receptor, TIR1 [27]. Therefore, miRNAs genome-wide pro ling during consistent heat stress in trees is necessary, so as to comprehensively understand the poplar mechanism in heat stress response and to examine the underlying useful genes in molecular breeding.
In this work, a photosynthesis regulatory network was established responding to heat stress by the transcriptome method coupled with weighted gene co-expression network analysis (WGCNA) analysis and graphic Gaussian model (GGM) algorithm. We show distinct transcriptional changes involving a highly interconnected hierarchically structured network and identi ed key nodes which response to heat stress in poplar. This study provides insights into the molecular mechanisms and systematic investigation of photosynthesis under heat stress and provides a new strategy for examining the genetic network and key nodes in plants.
Pro les of transcripts under heat stress in P. trichocarpa To investigate the transcriptional response under heat stress, we conducted high-throughput RNA-seq (Additional le 1; Average ~24 million reads). In total, 17,003 mRNA genes were detected to be highly expressed (FPKM > = 5, fold-change > 2 and P < 0.05), of which, 11,417 were differentially expressed at least at one of the six time points as compared to the control group (Table 1; Additional le 2). These DEGs (differentially expressed genes) are clustered into four co-expression modules by applying K-means clustering algorithm (Figure 2a). Genes of the four clusters exhibited distinct responsive patterns under the heat stress ( Figure 2b). For example, the transcript abundance of genes in cluster 4 increased at an early stage and reached a maximum level at 12 h, and then decreased after 36 h. The expression of genes in cluster 1 and 2 remained stable after stress (Figure 2b).
To examine the biological effect of heat responding genes, the functional enrichment analysis was carried out for identifying the over-represented gene ontology (GO) terms (Additional le 3). A total of 15 remarkably common GO terms were identi ed (P < 0.05), all of which were highly related to photosynthetic process (Supplementary Figure 1; Additional le 4). We noticed that a high percentage of genes from photosynthesis, cell wall, and chloroplast and cell division were enriched in the downregulated genes at 36 h. In total, 117 TF coding genes belonged to 26 families, which were speci cally expressed during these six treatment times (Additional le 5). We then used MapMan analysis, which groups heat-responsive genes into hierarchical functional categories based on their putative biological function. As a result, 168 genes were involved in photosynthesis regulation pathway, consisting of light reactions (92), Calvin cycle (41), and photorespiration (35). Of these 168 genes, 152 genes (90.48%) showed decreased expression patterns, as would be expected for decreasing the redox transient e ciency, assembling of Rubisco the transport and phosphorylation e ciency in photosynthesis. 76 from the 92 genes in light reaction that participated in key protein complexes were reduced, such as photosystem II (PSII),, photosystem I (PSI),, F-ATPase, and cytochrome b6f complex of redox chain (Cyt/b6f) (Figure 3a). In the photorespiration process, GLDC (Potri.001G144800), GLDP (Potri.006G229300), and other units taking part in converting glycine into methylene-H4-folate, together with SGAT (Potri.009G047700 and Potri.001G253300) in the peroxisome reaction with mitochondrion, were also decreased. Such nding indicated that glyoxylate metabolism had been repressed in the presence of heat stress, while enzymes had been deposited in mitochondria. We observed among genes downregulated in 36 h as well as an over-representation of genes involved in secondary metabolism (aspartate family amino acid, anthocyanin, organic acid, and phenylpropanoid biosynthesis) and primary metabolism (photosynthesis, carboxylic acid, tyrosine, amino sugar metabolisms and nucleic acid binding transcription factor). Above results suggested that heat leads to a reprogramming of primary and secondary metabolism in the poplar transcriptome putatively related to defense activation. Taken together, our result shows a dynamic transcriptional heat-responsive network, which had a negative effect on photosynthetic pathway.
Co-expression patterns responding heat stress in P. trichocarpa To investigate the trait-associated gene regulatory network (GRN) under heat stress, we constructed coexpressed gene matrix via weighted gene co-expression network analysis (WGCNA). A total of 17,003 genes were included in the analysis, and these genes exhibited a high expression dynamic across the

Heatresponsive genes and miRNAs involved in photosynthesis
To investigate the potential role of miRNAs in regulating photosynthetic, small RNA-seq (average ~13 million reads) were performed (Additional le 6). Our analyses revealed that the DEGs (differentially expressed genes) could be targeted by 163 miRNAs (Additional le 7). After investigating the regulation patterns, 54.63% miRNA/target pairs exhibited signi cant negative correlations (r = -0.73, P < 0.01). The highest numbers (8,718 genes) of DEGs were detected at 36 h ( Figure 2c, Table 1). Based on this criterion, we identi ed a total of 168 photosynthesis-related genes could be regulated by 159 miRNAs, of which, 88 are involved in light reaction and could be targeted by 72 miRNAs. Consistently, 50% miRNA/targets showed signi cant high negative correlations, suggesting that miRNAs play important role in regulating photosynthesis. In addition, in the Calvin cycle, 41 genes could be targeted by 33 miRNAs, and 41.3% showed negative correlations. For example, RBCS1A (Potri.015G062100) and RBCS3B (Potri.002G007100), protein subunits of Rubisco, were downregulated, while miR169 were up-regulated under heat stress ( Figure S7). These ndings suggest that a substantial number of the photosynthetic genes can be regulated by miRNAs under heat stress.

Genetic networks involved in photosynthesis response to heat stress
Combined with the weighted gene co-expression network (WGCNA) and MapMan analyses, we detected 77 differentially expressed genes (DEGs) involved in regulating photosynthetic pathways like electric transduction, carbon assimilation, and light reaction (Additional le 8). We identi ed 42 cis-regulatory motifs including some light-responsive elements that were enriched in the promoters of the 77 DE photosynthetic genes (P < 0.05), raising the possibility that some key regulators involved in this process. To further explore the key factors involved in regulating photosynthesis, we constructed a multilayers hierarchical gene regulatory networks that govern photosynthetic pathways. As a result, 77, 40 and 20, genes/miRNAs pairs were found to be located at the rst (bottom), second (middle), and third layers (top), respectively (importance score > 8.5, P < 0.05; Figure 6a). Consistently, 30 correlations are detected by a gene-gene correlation matrix based on the RT-qPCR result with a 0.72 (P < 0.001) correlation coe cient (Supplementary Figure 8).
Thirteen key TFs, such as Dof3, MYB22, ZF-HD3, ERF6, and MIKC-MAD4, were detected to regulate 74 downstream photosynthetic genes. HsfA1, a master regulator of transcriptional regulation under heat stress, induce the expression of the downstream transcription factors, such as DREB2A, MIKC-MAD, and NAC. In addition, DREB2A could induce the transcription of 30 downstream target genes. MiRNAs are also emerging as factors involved in regulating photosynthesis. Of these, miR1444d and miR396a could regulate 32 and 34 downstream photosynthetic genes, respectively, such as PETC, FD, PGRL1, and PETE, which are involved in electron transfer. Notably, the targets of some miRNAs were TFs, which ampli ed the effect of miRNAs through hierarchical regulation. For example, miR396 could regulate a MYB transcription factor (MYB22) and a bZIP DNA binding protein (bZIP28). Consistent with this, miR1444 were functionally signi cant miRNAs targeting polyphenol oxidase (PPO) genes for cleavage which started 2 nt upstream of precursor, resulting in differential regulation of PPO targets [34]. Notably, 19 downstream DEGs could be regulated by both key TFs and miRNAs. For example, MYB22 and miR396a were found to regulate 19 common DEGs (Figure 6b). Overall, the heat-responsive pathway in regulating photosynthesis is a multi-layered complex network which is co-controlled by TFs and miRNAs.

Discussion
Negative transcriptional effects of heat stress on photosynthesis Under global climatic change, elevated temperature is seen as the most serious threat to crop production [35][36][37]. Heat stress, whether transient or continuous, could cause biochemical, physiological, and molecular changes that adversely affect tree growth and productivity by reducing photosynthesis [38][39][40]. Photosynthesis in plants is composed of two interconnected biological functions, CO 2 transport and biochemical processes [41]. Studies on the ne temporal scale (i.e., hours) of the poplar response to heat are virtually rare, speci cally those dealing with understanding of changes in photosynthesis under heat stress. Here, we provide the rst glimpse into the regulation of photosynthesis under heat stress.
Previous studies have reported that heat stress induced cleavage and aggregation of PSII reaction center (RC) proteins, resulting in negative effects on photosynthesis [42,43]. High temperature could also inhibited photosynthesis by several mechanisms including deactivation of Rubisco, activation of NADPmalate dehydrogenase and an acceptor-side limitation imposed on photosynthetic electron transport in cotton, grapevine, tobacco and other plants [44][45][46]. Our study demonstrated that three important photosynthetic characteristics including transpiration rate (Trmmol), leaf internal CO 2 concentration (Ci), and stomatal conductance (Cond) were signi cantly decreased, suggesting that heat stress has permanent adverse effects on photosynthesis; however, no appreciable effect was observed for net photosynthesis (Pn). Previous studies investigated the main cause of reduced Pn and concluded that it may be caused by the decrease of Ci and Cond [11]. If both Cond and Ci simultaneously increase, increasing stomatal conductance will mainly limit Pn [40,42]. High temperature stress causes multi-step injuries to the photosynthetic machinery (PM), and the enzymes of the Calvin-Benson cycle are heat labile [42]. Therefore, the carbon absorption system shows high sensitivity to the increased temperature, which is greatly suppressed under moderate heat stress [12,13]. The reduced Rubisco activity in the presence of moderate heat stress is related to photosynthetic loss [43]. Heat will induce alterations in Chl-protein complex structures and inactive enzyme activities, which have constituted the heat stress effects, eventually contributing to decreasing Pn. Heat stress at 35-45 ℃ can result in generally reorganized thylakoid membranes and reduced membrane stacking [47]. The above structural changes resulted in a modest decrease of Cond and Ci, which might cause Pn to decrease.

A model for the heat stress-responsive pathway in regulating photosynthesis
Previous studies have demonstrated extensive changes in the transcriptome under heat stress [1,3,42]. Nonetheless, the alterations of gene regulation network in woody plant photosynthetic tissues have been rarely mentioned. Here, by using dynamic transcriptome data, this paper had demonstrated numerous negative and positive transcriptional regulators, which might function as the transcriptional cascade to activate downstream target gene expression in the photosynthetic pathways: light reaction, photorespiration, and Calvin cycle. Indeed, several cis-elements were co-occurred in the promoter regions of the down-stream photosynthetic genes. A total of 42 cis-elements were enriched in the promoters of 77 downstream photosynthetic DEGs, such as the TATA and CAAT boxes, stress responding elements (BIHD1OS and EBOXBNNAPA), light responding elements (GATABOX and SORLIP1AT), hormone responding elements (WRKY71OS and MYB1AT), organ speci c elements (POLLEN1LELAT52 and GGNCC), the phosphate-starvation responding element, and the sulfur responding element. These results thus provide potential for studying the heat-responsive genes involved in photosynthesis.
The high temporal-resolution that transcriptome pro ling data generated in our study provided good opportunity to identifying key regulators, which is also highly informative for inferring gene function and understanding the genetic control of regulating photosynthesis under heat stress. In the present study, we identi ed 4 TFs and 6 miRNAs, which are key players in regulating photosynthesis. In upstream regulating cascade, heat stress can derepress the signal transduction pathway through promoting proteolysis of the negatively regulated miRNAs through interaction for releasing TFs [3]. MYB22, a key TF in the regulation network, known to be involved in the maturation program in development under heat stress [48,49]. MiR396, also a key node detected in our network, could repress the expression of PIF1 and further impact the transcription of SBPase [50,51]. For downstream target genes, LHCA, LHCB, RBCS, and ATPase were repressed under heat stress across different treatment times, supporting the idea that upstream receptors and signaling components may respond at different heat stress treatment times under. Thus, taking advantage of the high temporal-resolution transcriptome, we could propose a model for the heat stress-responsive pathway in regulating photosynthesis (Figure 7).

A novel approach for gene regulatory network construction
It has been proved that, it is useful to apply gene expression pro ling data in constructing co-expression network, and in performing network decomposition and analysis in biological study. Combining the expression-traits association method, cis-elements enrichment, and multilayered hierarchical gene regulatory network (ML-hGRN) construction approaches to discovery key players controlling photosynthetic genes is particularly important. Weighted correlation network analysis (WGCNA) is distinctly different from previously mentioned methods in that the former de nes a network persistently linking all variables before clustering variables with high co-expression within the modules de ned exibly [52][53][54]. It could explicitly de ne all of the relationships between metabolites at multiple scales, facilitates identifying essential downstream processes that are tightly related to traits. In our study, we rst identi ed 1,548 downstream heat-responsive genes. Then, based ciselement, we further reduced false-positives. Following stringent ltration, many genes are excluded from the eventual datasets; which do not stand for the direct active players in this pathway. Thus, 77 downstream photosynthetic genes and 442 upstream regulators, including 279 TFs and 163 DE miRNAs remained. Finally, based on the downstream target genes, and upstream regulators (TFs and miRNAs), a multilayered hierarchical gene regulatory network (ML-hGRN) was constructed. In this ML-hGRN, MYB22, miR396a, miR1444d, DREB2A, and ZF-HD3 were identi ed as key nodes, which were supported by previous studies. For example, DREB2A were reported to bind to a cis-acting dehydration-responsive element (DRE) and activated the expression of Hsp70, AREB1, and ATEM6 in response to drought and salt stress [24,55]. MiR396a, regulators of stress resistant processes, implicated to be key players in response to heat stress [33,50]. Therefore, this multifaceted approach identi ed key genes in regulating photosynthesis, and serves as a model for interrogating complex signaling networks.

Conclusions
Heat stress made extensive changes in the transcriptome and caused biochemical, physiological, and molecular changes that adversely affect plant growth. Although a large amounts of genes have been identi ed in many plant species, little attention has been paid to changes in gene regulation network in photosynthetic tissues of woody plants. Here we combined the expression-traits association method, ciselements enrichment, and multilayered hierarchical gene regulatory network (ML-hGRN) approaches to discovery key players controlling photosynthetic genes. A total of 10 key nodes (4 TFs, 6miRNAs) were identi ed which regulating 77 essential photosynthetic pathway genes under heat stress. Moreover, the expression pattern of 77 photosynthetic genes at different treatment points in response to heat stress revealed that different members might play speci c roles in responding heat stress. Interestingly, the expression pro ling of photosynthetic genes showed similar down-regulating pattern under heat stress across different time points, suggesting stress-induced genes might have multiple functions in downregulating photosynthesis and responding heat stress. Furthermore, the heat-responsive pathway in regulating photosynthesis is a multi-layered complex network which is co-controlled by TFs and miRNAs. Importantly, MYB22, miR1444 and miR396a might make signi cant roles in responding to abiotic stresses, particularly in heat and drought. Taken together, the results imply a functional role for these key genes in mediating photosynthesis responding to abiotic stress, demonstrating time-course transcriptome-based regulatory network construction will facilitate further the genetic network and key nodes examining in plants.

Plant Materials
Ramets of P. trichocarpa (one genotype) were used and collected by Dr. Yiyang Zhao, from a greenhouse located in Beijing Forestry University, Beijing, China (40°0′N, 116°08′E). The specimen was deposited in the Beijing Forestry University Herbarium (voucher 00080269). The samples were treated according to the method of Zhang et al. (2010). Ramets were transplanted to the nutrition blocks supplemented with common garden soil within a period of 30 days. Subsequently, samples that had identical growth degree were transplanted to the nutrition pots (with the diameter of 30 cm) supplemented with identical amount of common garden soil separately. Then, P. trichrocarpa ramets at the age of 1 year had been planted at the 16 h light and 8 h dark cycle. Subsequently, seedlings had been exposed to 40 °C temperature for with 4,8,12,24,36, and 48 h, and samples that grew under 25 °C had been utilized as controls. The relative humidity was determined at 50% ± 1, which was maintained in the measuring process. All treatment groups, such as control, included 3 biological replications (ramets).

Measurement of photosynthetic rates
From August 18th-24th, 2017, the No.4 completely expanding leaves were collected from all the three ramets under every treatment to measure the photosynthetic rate through the portable photosynthesis system (LI-6400; Li-Cor Inc., Lincoln, NE, USA). In order to fully induce photosynthesis, every sample was subjected to 30 min of illumination at saturated photosynthetic photon ux density (PPFD) from the light emitting diode (LED) light source prior to measurements. Afterwards, the transpiration rate (Trmmol), net photosynthetic rate (Pn), stomatal conductance (Cond), and intercellular CO2 concentration (Ci) in leaf had been determined simultaneously. Each measuring parameter followed the methods mentioned by Chen et al. [56].
Construction, processing, and annotation of the P. Following analysis, those samples that had the OD 260/280 ratio of 1.9-2.1 and the OD 260/230 ratio of 1.8-2.2 would be selected for subsequent use. Paired-end sequencing was performed by Novogene (Beijing, China) on an Illumina HiSeq 4000 platform (Illumina), generating 150-bp paired-end reads using TruSeqPEClusterKitv3-cBot-HS (Illumina) according to the manufacturer's instructions and the strandspeci c libraries were sequenced. Quality control and other downstream analyses were performed on the cleaned data. In this step, clean data (clean reads) had been acquired through eliminating adapter-and ploy-N-containing reads, as well as the low-quality reads out of the raw material. The clean reads had then be mapped into reference genome of P. trichrocarpa with Bowtie2 program. The ltered high-quality reads were mapped on the P.trichrocarpa genome using TopHat (v2.0.0) with default parameters. Feature Countsv1.5.0-p3 had been utilized for counting reads that were mapped into every gene; afterwards, FPKM was computed for all genes according to gene length and reads that were mapped into the speci c gene. The mapped results were then processed using Cu inks (v2.0.2) for obtaining FPKM for every Populus gene of every sample with default parameters. Then, the correlations among biological replicates had been examined through computing the Spearman correlation coe cient. Cuffdiff was used to determine the differential expression among various samples. Genes that had > = 2-fold change and the P-value of < = 0.05 [adjusted with the false discovery rate (Q-value)] had been deemed as the signi cantly differentially expressed genes (DEGs).
Small RNA (sRNA) library construction, annotation and processing for P. trichrocarpa Total RNA sample, which was identical to that utilized in transcriptomic sequencing, was used for constructing sRNA libraries. In brief, the sRNA pieces (18-30 nt) had been isolated, which was then subjected to polyacrylamide gel electrophoresis (PAGE) for puri cation. Afterwards, the sRNA segments had been ligated to 3' and 5' adaptors using T4 RNA ligase. Then, the adaptor-ligated sRNAs had been subjected to transcription into single strand cDNA, followed by PCR ampli cation for 12 cycles. Later, samples coded by index were clustered using the cBot Cluster Generation System by the TruSeq SR Cluster Kit v3-cBot-HS (Illumia) in accordance with manufacturer protocols. Following cluster production, all library preparations had been sequenced through the Illumina Hiseq 2500 platform to generate the 50 bp single end reads.
Identifying differentially expressed miRNAs and analyzing the speci c targets In each sample, each read was standardized to be the transcripts per million (TPM) based on all sequenced reads. Typically, miRNAs that had the log 2 fold change of > = 1.0 or the log 2 fold change of < = −1.0, FDR < = 0.05, and P < 0.05 had been deemed as the signi cantly differentially expressed miRNAs. Speci cally, fold change was calculated as follow, fold change = |TPM ratio| (treatment/control). Then, psRNATarget (http://plantgrn.noble.org/psRNATarget/) with accessibility of parameter target at 25.0, together with the complementarity scoring length at 20, had been utilized for predicting the candidate target genes for a total of 179 miRNAs. Typically, targets were identi ed using the transcript library of P. trichocarpa preserved within Phytozome database v.3.0 (http://www.phytozome.net/poplar). Moreover, the scoring systems for the mismatches during the interaction between miRNA and mRNA had been assigned below: 0 suggested a complementary pair; 0.5 indicated G:U wobble; 1 stood for other mismatches; 2 represented every opening gap; while 0.5 was indicative of the extending gap [57]. Moreover, the miRNA/mRNA duplexes that had the score of < = 2.0 had been regarded to be the highly con dent target-miRNA interactions. Accordingly, candidate targets that had the fold change of > = 2, Q < 0.05 and P < 0.01 had been deemed with signi cant differential expression. Additionally, to examine the relationship of expression between miRNAs and target genes, pairwise Pearson correlations between targets and miRNAs had been computed to explore the association of miRNA expression with target gene expression. All P-values had been adjusted by the false discovery rate (FDR; P < 0.01, Q < 0.05). The R software (R Development Core Team, 2012) was used for calculation.

GO and pathway enrichment analysis
The AgriGO analysis approaches (http://bioinfo.cau.edu.cn/agriGO/) were utilized for enrichment analysis of gene set. The precise annotation data were obtained based on the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg), as well as PopGenie (http://popgenie.org).
Then, the enrichment P-value was computed for every represented GO term, which was also adjusted according to the Benjamini Hoschberg error correction approach. Afterwards, GO terms that had adjusted (after FDR adjustment) P≤0.05 had been regarded as signi cantly enriched. Furthermore, the MapMan (Version: 3.5.1) categories were used for pathway enrichment analysis of various gene sets (signi cance value < = 0.05).

Co-expression network analysis for construction of modules
For co-expression network analysis, the weighted gene co-expression network analysis (WGCNA) [58] package had been employed. According to the log 2 (1 + FPKM) numbers, the pairwise Spearman correlation coe cient matrix had been generated among different gene pairs, which was then converted into the adjacency matrix (the connection strength matrix) through the following formula: connection strength (namely, adjacency value) = |(1 + correlation)/2| β , where β represented the correlation matrix soft threshold and it rendered a higher weight to the strongest correlation in the meantime of keeping connectivity between genes. Additionally, the β-value at 12 had been chosen according to scale-free topological criteria from Horvath and Zhang (2005). Afterwards, the obtained adjacency matrix would be transformed into the topological overlap (TO) matrix (TOM) using the TOM similarity algorithm; then, genes had been clustered hierarchically according to the TO similarity. Hierarchal clustering dendrogram was cut by dynamic tree-cutting algorithm, and then modules would be de ned following branch decomposition/combination, so as to obtain the xed cluster quantity [58]. In every module, PCA was computed through the summary pro le (namely, module eigengene, ME). Besides, modules that had greater TO values (average TO of all genes within one speci c module) than those of modules consisting of genes selected randomly would be preserved. Later, Cytoscape [59,60] was used for GO enrichment analysis for all modules, as described above.
Identi cation and co-expression analysis of cisregulatory elements and TFs The transcription modules were predicted in accordance with the methods by Belmonte et al. (2013) after certain modi cations. In order to examine over-represented DNA motifs within the heat responding photosynthetic DEGs, promoter sequences (2 kb) in front of the 5 terminal in every annotated transcriptional unit during Populus annotation had been collected. Then, the identi ed promoter elements of plants, together with the corresponding annotations, had been derived from PlantPAN 2.0 [61].
Moreover, the promoter sequences from each Populus gene had been utilized to be the background, and hypergeometric test was carried out to detect the over-represented motifs. A co-expression analysis (Pearson correlation) was then performed using the Hmisc 4.1.0 package in R for detecting the correlation networks of gene expression, and all criteria below were used for gene ltering: |r| > = 0.6 and P < = 0.05. The networks were drawn using Cytoscape [59,60].
Multilayered hierarchical gene regulatory network (ML-hGRN) construction using a backward elimination random forest (BWERF) algorithm The sequences of the 77 genes (Additional le 8) known to be directly involved in photosynthesis which were classi ed as light reaction, Calvin cycle and photorespiration based on MapMan and WGCNA analyses, and the PlantPAN2.0 was utilized to evaluate the cis-regulatory elements in every promoter sequence (2 kb in length) [61]. O the basis of 80% con dence, the above motifs had been utilized in identifying the 279 TFs that potentially targeted the photosynthetic genes using PlantPAN2.0.
Transcriptional pro les of these TFs and photosynthesis genes were obtained from the transcriptome data, and those genes were retained for further analysis. A ML-hGRN had been established according to the BWERF algorithm, among which, 77 photosynthetic genes were in bottom layer, while 163 miRNAs and 279 TFs were within regulating layer [62]. Additionally, a three-layered GRN had also been constructed using 32 miRNAs and 18 TFs potentially regulating the 77 photosynthesis genes either directly or indirectly. Cytoscape was used to draw the network [59].

Quantitative realtime PCR (qRTPCR) analysis
The 7500 Fast Real-time PCR System (Applied Biosystems, Waltham, MA, USA) was used for RT-PCRs by Light Cycler-FastStart DNA master SYBR Green I kit (Roche). The Primer Express 3.0 (Applied Biosystems, Waltham, MA, USA) was employed to design each primer pair of purative gene (Additional le 9). Each RT-qPCR ampli cation was performed according to the standard reaction procedure in triplicate biological and technical replications. Then, the ampli ed fragment speci city had been examined by the produced melting curve. The Opticon Monitor Analysis Software 3.1 tool (Bio-Rad) was adopted to analyze the produced real-time data. To measure the expression of photosynthetic genes, the poplar Actin gene (accession no. EF145577) had been utilized for internal reference. The relative miRNA level would be measured and normalized based on 5.8S rRNA according to the method from Song et al. (2013). Moreover, the total RNA was reversely transcribed into the cDNA templates used in reactions.