Comparative transcriptomics and bioinformatics analysis of genes related to photosynthesis in Eucalyptus camaldulensis

Plant material

In early April 2020, leaf samples of E. camaldulensis were collected from Southern China Experiment Nursery in Zhanjiang, Guangdong, China (21°15′30.69″N, 110°06′41.95″E). The young leaves (1-month-old) and mature leaves (4-month-old) of three healthy E. camaldulensis plants in the same direction and same height of second canopy were collected from the same half-sib families located on the Morehead R. Queensland (15°15′S, 143°34′E) with three biological replicates. The samples were kept in liquid nitrogen for RNA extraction.

Determination of photosynthetic indexes

To measure the photosynthetic indexes of E. camaldulensis leaves, a portable photosynthetic apparatus (LI-6400XT) was used under natural light (from 10:00 to 12:00), which included the net photosynthetic rate (Pn), cond (Gs), trmmol (Tr), and leaf water use efficiency (WUE) = Pn/Tr. The stability of the environment was maintained by setting the light intensity at 1200 µmol m⁻² s⁻¹, the flow rate at 500 µmol s⁻¹, and the CO₂ concentration at the current atmospheric CO₂ concentration ((400 ± 20) µmol mol⁻¹), and the leaf temperature at 25 ± 1 °C. A total of three branches were selected, and three young and mature leaves were measured respectively in different parts of each branch with three biological replicates.

Measurement of the light response curve

LI-6400XT portable photosynthetic apparatus was used to measure the light response curve of E. camaldulensis leaves from 10:00 to 16:00. The experimental apparatus was set as an open-air circuit keeping CO₂ at 400 µmol mol⁻¹, and the gas flow rate was 500 µmol s⁻¹, and other parameters were set as default. A total of 17 photosynthetically active radiation (PAR, µmol m⁻² s⁻¹) were set using LED red and blue light sources, which included 2000, 1800, 1600, 1400, 1200, 1000, 800, 600, 400, 200, 150, 100, 80, 60, 40, 20 and 0 µmol m⁻² s⁻¹. The net photosynthetic rates of leaves were measured at different light intensities. Before determining the photosynthetic rate, they were induced for 2 min under the light intensity of 1200 µmol m⁻² s⁻¹, and the waiting time for each light intensity point was 120–200 s. A rectangular hyperbola model was used to draw the light response curve (Baly, 1935). ${\displaystyle A_{\mathrm{n}}\left(I\right)=aI{A}_{\max }/\left(aI+A_{\max }\right)-R_{\mathrm{d}}.}$ A_n(I) represents the net photosynthetic rate, I represents the light intensity, a represents the initial slope of the light response curve, A_max represents the maximum net photosynthetic rate, R_d represents the dark respiration rate. The analysis methods of one-way ANOVA and Duncan’s were used in SPSS19.0.

RNA extraction and transcriptome sequencing

Total RNA was extracted using FastPure Universal Plant Total RNA Isolation Kit (Vazyme, Nanjing, China). The quantification of RNA was performed using Nanodrop™ 1000 (Thermo Fisher Scientific, Waltham, MA, USA). Agilent 2100 Bioanalyzer was used to measure the 28S/18S and RIN values. The quality and integrity of the RNA samples were estimated using agarose gel electrophoresis, and samples were resolved on 1% agarose gel. The enriched mRNA was reversely transcribed to form double-stranded cDNA, and the double ends of cDNA were repaired, and the adaptor was added for PCR amplification to construct the on-machine library. RNA libraries were constructed using the Illumina HiSeqTM 4000 platform at Gidio company. The statistical power of this experimental design, calculated in RNAseqPS is 0.446 (Guo et al., 2014).

Sequencing data processing and analysis

TopHat (v2.0.10) was used to compare the reads with E. grandis reference genome (NCBI gcf_000612305.1). The assembly of sequences was performed by Cufflinks using reference annotation-based transcripts (RABT). The obtained transcript sequences were compared using BLAST with SwissProt, Gene Ontology, Cluster of Orthologous Groups of proteins, and Kyoto Encyclopedia of Genes and Genomes databases to annotate the function of the transcript.

The differentially expressed proteins were mapped to each term of the GO databases, and calculated the number of proteins in each term, and obtained the list of proteins with a certain GO function and the number of proteins. Hypergeometric tests were used to identify GO entries that were significantly enriched in differentially expressed proteins compared with the whole background protein, and the p-value of this hypothesis test could be calculated as follows: ${\displaystyle P=1-\sum _{i=0}^{m-1}\frac{\left(\begin{array}{c}\hfill M\hfill \\ \hfill i\hfill \\ \hfill \hfill \end{array}\right)\left(\begin{array}{c}\hfill N-M\hfill \\ \hfill n-i\hfill \\ \hfill \hfill \end{array}\right)}{\left(\begin{array}{c}\hfill N\hfill \\ \hfill n\hfill \\ \hfill \hfill \end{array}\right)}.}$

N was the number of proteins with GO annotation among all background proteins; n was the number of differentially expressed proteins in N; M was the number of proteins annotated as a specific GO term among all background proteins; m is the number of differentially expressed proteins annotated as a specific GO term. After the calculated p-value was corrected by Benjamin&Hochberg(BH method), the threshold was taken as correct-p-value ≤0.05. The major biological functions of differentially expressed proteins could be determined by GO functional significance enrichment analysis.

Pathway significant enrichment analysis took KEGG pathway as the unit, and hypergeometric test was used to find out the pathways that were significantly enriched in differentially expressed proteins compared with background proteins. The p-value of this hypothesis test could be calculated as follows: ${\displaystyle P=1-\sum _{i=0}^{m-1}\frac{\left(\begin{array}{c}\hfill M\hfill \\ \hfill i\hfill \\ \hfill \hfill \end{array}\right)\left(\begin{array}{c}\hfill N-M\hfill \\ \hfill n-i\hfill \\ \hfill \hfill \end{array}\right)}{\left(\begin{array}{c}\hfill N\hfill \\ \hfill n\hfill \\ \hfill \hfill \end{array}\right)}.}$

N was the number of proteins with pathway annotation among all background proteins; n was the number of differentially expressed proteins in N; M was the number of proteins annotated as a specific pathway among all background proteins; m was the number of differentially expressed proteins annotated as a specific pathway. After correction for multiple tests, pathways with q-value ≤0.05 were defined as those significantly enriched in differentially expressed proteins. Here, q-value is the p-value after FDR correction. The significant enrichment of pathway can identify the main biochemical metabolic pathways and signal transduction pathways involved in differentially expressed proteins.

Cuffdiff software (Trapnell et al., 2012) and EdgeR (Robinson, McCarthy & Smyth, 2010) were used to analyze the transcriptional expression and differences between mRNA sample pairs, respectively. DEGs were selected with a p-value <0.05 and —log2 fold change (log2FC)—>1 for subsequent analyses.

Gene analysis

DEGs related to photosynthesis were screened out based on the data of and transcriptome. To analyze the relative molecular weight, theoretical isoelectric point, instability coefficient, hydrophilicity index as well as lipid solubility index of the gene, the ProtParam in ExPASy (http://web.expasy.org/protparam/) were employed. The secondary and tertiary structures of protein were analyzed using the SOPMA (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html) and SWISS-MODEL (https://swissmodel.expasy.org/interactive#seque-nce), respectively.

RNA reverse transcription and QRT-PCR

M-MLV reverse transcriptase was adopted for synthesizing 2 µg RNA in line with Evo M-MLV RT Kit instructions (Accurate Biotechnology, China), and followed by a 1:10 dilution for subsequent experiments. The photosynthesis-related genes of E. camaldulensis were identified as verification genes. The total qRT-PCR reaction comprised of 20 µL containing 2 ×SYBR Green Pro Taq HS Premix (Tli RNaseH Plus) 10 µL, ROX Reference Dye 0.4 µL, forward primer (10 µmol L⁻¹) 0.8 µL, reverse primer (10 µmol L⁻¹) 0.8 µL, cDNA 2 µL, and RNase Free ddH₂O 6 µL. The thermal profile was comprised of two segments: 30-S at 95 °C; 5-S denaturation at 95 °C, and 30-S annealing at 60 °C for altogether 40 cycles. Each reaction was repeated three times, and the actin gene was selected as a reference. Primer Express 2.0 Software (PE Applied Biosystems, USA) was applied in primer designing with default parameters. Table S1 displays the sequences of all the primers. The 2^−ΔΔCt method (Livak & Schmittgen, 2001) and Microsoft Excel software were used to analyze the obtained data.

Photosynthesis analysis of the mature and young leaves in E. camaldulensis

Photosynthesis analysis revealed that mature leaves’ photosynthetic rate, trmmol, cond, and leaf water use efficiency were significantly higher than those of young leaves at a light intensity of 1200 µmol m⁻² s⁻¹ (Fig. 1). A photosynthesis rate is also called as a photosynthetic intensity, which is a measure of the intensity of photosynthesis. The photosynthetic rate can be expressed as carbon dioxide absorbed or oxygen released per unit time and per unit leaf area. The light response curve reflects the variation of plant photosynthetic rate with increasing light intensity; the net photosynthetic rate increased rapidly from 0 to 600 µmol m⁻² s⁻¹. However, the net photosynthetic rate decreased slower at light intensities of 600–1600 µmol m⁻² s⁻¹. A higher net photosynthetic rate was observed for mature leaves than for young leaves (Fig. 2).

Differentially expressed genes related to photosynthesis in mature and young leaves of E. camaldulensis

We carried out a transcriptomic analysis on both young and mature leaf samples of E. camaldulensis to explore potential target genes. The overall raw/clean reads within each sample ranged between 6,914,804,063 and 8,686,921,500. This work also aligned sequence reads into the reference genome of E. grandis, and results indicated that >81% were map-able. The >92% Q30% and 48% GC concentrations suggested high-quality transcriptomic results for the subsequent analyses (Table 1).

The results underlined that there were 18,443 differentially expressed genes (DEGs) in young and mature leaves in E. camaldulensis. A GO enrichment analysis indicated that DEGs were mainly enriched in the following processes: single-multicellular organism process, multicellular organism process, anatomical structure development, single-organism developmental process, multicellular organism development, and developmental process (Fig. 3). Furthermore, KEGG enrichment indicated that DEGs were mainly enriched in plant hormone signal transduction, carbon metabolism, photosynthesis, and carbon fixation in photosynthetic organisms (Fig. 4). We focused on the DEGs relating to photosynthesis, as mature leaves have a higher net photosynthetic rate than young leaves (Figs. 1 and 2).

A total of 18 important DEGs were identified in E. camaldulensis related to photosynthesis at different periods. Among ten ribose-5-phosphate isomerases genes (RPI), six were up-regulated in the mature leaves compared to young leaves, while the other four were down-regulated. There were four differentially expressed rubisco methyltransferase (RBCMT); two of them were up-regulated as well as two were down-regulated in the mature leaves compared to the young leaves. Furthermore, two ribulose-phosphate-3-epimerase genes (RPE) were found, and both of these were up-regulated in the mature leaves in comparison to the young leaves. On the other hand, two triosephosphate isomerase genes (TIM) were down-regulated in the mature leaves relative to the young leaves (Table 2).

Gene analysis bioinformatics analysis of genes related to photosynthesis

In order to further analyze the 18 DEGs related to photosynthesis in E. camaldulensis, bioinformatics methods were used to analyze them. Among the RPI genes of E. camaldulensis, RPI2 was the longest amino acid encoded protein with a maximum relative molecular mass of 29,632.15 D, whereas the minimum relative molecular mass was RPI9 (3852.50 D). RPI2, RPI3, and RPI5 genes encoded for basic amino acids, whereas others encoded for acidic amino acids. A minimum instability index of 0.46 was observed in the case of RPI3. According to the hydrophilic index, the proteins encoded by RPI1, RPI5, RPI8, and RPI9 genes were hydrophobic in nature, whereas others were hydrophilic proteins. The aliphatic index of RPI was found between 83.33 and 108.11. In the RBCMT genes of E. camaldulensis, RBCMT2 was the longest amino acid encoded protein with a maximum relative molecular mass of 54957.68 D, while the minimum relative molecular mass was RBCMT4 (9489.89 D). While RBCMT1 encoded basic amino acids, others encoded acidic amino acids. The minimum instability index was 46.48 (RBCMT3). According to the hydrophilic index, the proteins encoded by RBCMT genes were hydrophobic in nature. The aliphatic index of RBCMT was found between 86.75 and 114.10. RPE1 and RPE2 encoded for the basic and acidic amino acids, respectively. The protein encoded by the RPE1 gene was hydrophobic in nature, whereas the RPE2 gene encoded for hydrophilic protein. TIM1 and TIM2 encoded acidic amino acids as well as hydrophilic proteins (Table 3).

The protein structures of the RPI, RBCMT, RPE, and TIM gene consisted of α-helix, extension chain, random coil, and β-angle structure (Table 4 and Fig. 5). These results showed that the proportion of RPI protein α-helix structures was usually higher than others, with the lowest is β-angle structure. Similarly, the proportion of RBCMT protein α-helix and random coil structure was also usually higher than that of protein with other structures. The highest proportion of α-helix, as well as the lowest ratio of β-angle, was found in the secondary structure of RPE and TIM proteins.

Q RT-PCR analysis of key genes related to photosynthesis

Ten of the 18 target genes selected for analysis were over-expressed in mature leaves, namely RPI1, RPI2, RPI3, RPI4, RPI5, RPI6, RBCMT1, RBCMT2, RPE1, and RPE2. The higher expression of these genes in mature leaves than younger ones was verified by qRT-PCR (Fig. 6A). On the other hand, eight target genes were down-regulated in the mature leaves viz., RPI7, RPI8, RPI9, RPI10, RBCMT3, RBCMT4, TIM1 as well as TIM2 and the qRT-PCR verified that their higher expression in young leaves than mature leaves (Fig. 6B). A consistent pattern of results was found between the results of the qRT-PCR and the RNA sequencing.

Compare photosynthesis between young and mature leaves in E. camaldulensis

The rate of photosynthesis is closely related to the leaf age, and it changes with the age of the leaves. The photosynthetic rate showed a single peak curve from leaf initiation to senescence and wilting (Smith et al., 1997). In this study, the photosynthetic rate and trmmol of mature leaves were found to be higher than those of young leaves in E. camaldulensis. Similar findings were observed in Chinese fir, where the light absorption capacity and transformation capacity of old leaves were reported to be higher than those of young leaves (Li et al., 2017). This phenomenon may be possible due to the higher vulnerability of chlorophyll in young leaves to intense light. An increase in the net photosynthetic rate leads to higher stomatal conductance followed by more CO₂ into carboxylated parts of the plant cells for photosynthesis. The increase of stomatal conductance also accelerated the loss of water, leading to the continuous increase of transpiration rate (De Souza et al., 2020). The light response curve reflects the change in net photosynthetic rate with the change of light intensity (Poorter & Navas, 2003), which provides the estimate of photosynthetic rate, light saturation point, light compensation point as well as other important ecological and physiological parameters of plants (Cannell & Thornley, 1998). Therefore, in this study, we analyze the light response curves of young and mature leaves of E. camaldulensis to understand the photosynthetic efficiency and growth of eucalyptus.

Key genes associated with photosynthesis in E. camaldulensis

The high-throughput sequencing provided huge data with high quality, which was suitable for the transcriptome study of Eucalyptus species, and improved our understanding about the molecular mechanisms in Eucalyptus (Morozova & Marra, 2008; Kadam et al., 2013; Xiao et al., 2020). Although, the role of Rubisco and PEP (two key enzymes in photosynthesis dark reactions) have been widely explored, there is limited information available on ribose-5-phosphate isomerase (RPI), rubisco methyltransferase (RBCMT), ribulose-phosphate 3-epimerase (RPE), and triosephosphate isomerase (TIM) genes. Enzymes play an important role in the dark reaction of photosynthesis; therefore, RPI, RBCMT, RPE, and TIM genes were selected as the research objects.

RPI is a highly conserved protease that is ubiquitous in many organisms and plays a central role in the pentose phosphate pathway (PPP). It is involved in the reversible isomerization of ribose-5-phosphate (R5P) and ribulose 5-phosphate (Ru5P) between prokaryotes and eukaryotes as well as in the Calvin cycle of carbon dioxide fixation in plants (Wamelink et al., 2010; Gontero, Cárdenas & Ricard, 1988). According to the reports, RPI is an essential enzyme in the Calvin cycle in peas and spinach (Woodruff & Wolfenden, 1979; Park & Anderson, 1973). The photosynthetic rate of mature leaves was higher than that of young leaves in E. camaldulensis, genes such as RPI1, RPI2, RPI3, RPI4, RPI5, and RPI6 were up-regulated in the mature leaves. Among the RPI genes of E. camaldulensis, RPI2 was the longest amino acid encoded protein. RPI2, RPI3, and RPI5 genes encoded for basic amino acids. The proteins encoded by RPI1 and RPI5 genes were hydrophobic in nature. The results indicate that the high expression of RPI genes in mature leaves could produce more PRI enzymes for photosynthesis.

Rubisco is a bifunctional enzyme and is widely available in the chloroplast matrix of plants. It can catalyze both the carboxylation reaction of C₃ for photosynthesis and the oxygenation reaction of C₅ for photorespiration. In photosynthesis, Rubisco was present at the intersection of carbon oxidation and carbon reduction, and RBCMT was the rubisco methyltransferase (Miyazawa et al., 2020). The expression of RBCMT1 and RBCMT2 was positively correlated with the net photosynthetic rate. RBCMT2 was the longest amino acid encoded protein in E. camaldulensis. RBCMT1 encoded basic amino acids. According to the hydrophilic index, the proteins encoded by RBCMT genes were hydrophobic in nature.

RPE can play a vital role in the development of the NADPH pool, and PPP, which can convert monosaccharides, such as glucose, into the nucleotide precursor of pentose sugars. Additionally, RPE can convert the ribulose-5-phosphate into xylulose 5-phosphate in the Calvin cycle (Lykins et al., 2018). The photosynthetic rate of mature leaves was higher than that of young leaves in E. camaldulensis. RPE1 and RPE2 were up-regulated in mature leaves, which indicates that these may be the key genes to enhancing photosynthesis. RPE1 encoded for the basic amino acids, and RPE2 encoded for the acidic amino acids. The protein encoded by the RPE1 gene was hydrophobic in nature, whereas the RPE2 gene encoded for hydrophilic protein. Furthermore, RP7, RP8, RP9, RP10, RBCMT3, RBCMT4, TIM1, and TIM2 genes were highly expressed in young leaves of E. camaldulensis, indicating that they may be involved in the negative regulation (Katrin et al., 1994). These results provide a basic structure for understanding the RPI, RBCMT, RPE, and TIM genes and exploring the functions of genes for further research in E. camaldulensis. Through genetic transformation, enhancing photosynthetic efficiency or productivity by using these genes to ensure the better afforestation quality of eucalyptus plantations and to generate greater social, ecological, and economic value.