Amino acids profiling and transcriptomic data integration demonstrates the dynamic regulation of amino acids synthesis in the leaves of Cyclocarya paliurus

Plant materials and sampling

C. paliurus leaves at three different developmental stages were collected on 01 May, 2018, in Zhuzhang village, Longquan city, Zhejiang Province, China (E118°48′28″, N28°5′57″). The S1 stage was during the smallest expanded leaf without curly leaf, while the S3 stage was during the largest expanded leaf with full leaf enlargement and full leaf thickness; and S2 stage was an intermediate developmental stage between S1 and S3 stages (Fig. 1). C. paliurus leaves were sampled separately on the same tree at the same time, then leaves were divided into three developmental stages according to the leaves’ developmental status and size (Lin et al., 2020). Leaves at S3 stage and S2 stage were significantly larger and longer than those at S2 stage and S1 stage, respectively (Table S1). To obtain sufficient leaves, three samples (three plants) were randomly selected and mixed as one replicate. Three biological replicates (nine plants) were used in this study. The samples were frozen immediately in liquid nitrogen and stored at −80 °C for subsequent metabolite extraction and transcriptome sequencing.

Metabolite extraction and profiling

Metabolite extraction and identification were conducted by Wuhan Metware Biotechnology Co., Ltd. (Wuhan, China) as previously described in Lin et al. (2020). The freeze-dried leaf was crushed using an MM 400 mixer mill (Retsch, Baltimore, MD, USA) with a zirconia bead for 1.5 min at 30 Hz. One hundred milligrams of C. paliurus leaf powder was extracted with 1 mL 70% aqueous methanol solution overnight at 4 °C. After centrifugation at 10,000g for 10 min, the extracts were absorbed (CNWBOND Carbon-GCB SPE Cartridge, 250 mg, 3 mL; ANPEL, Shanghai, China) and filtered (SCAA-104, 0.22 μm pore size; ANPEL, Shanghai, China) (Wang et al., 2019). All the samples were mixed and used as the quality control sample (QC) to monitor the precision of the assay.

Next, the sample extracts were analyzed using an LC-ESI-MS/MS system (UPLC, Shim-pack UFLC SHIMADZU CBM30A system; MS, Applied Biosystems 6500 triple quadrupole-linear ion trap mass spectrometer (Q TRAP)) with column (Waters ACQUITY UPLC HSS T3 C18, 1.8 µm, 2.1 mm × 100 mm), equipped with an ESI Turbo Ion-Spray interface, operating in both positive and negative ion mode and controlled using Analyst 1.6 software (AB Sciex, Redwood, CA, USA). Solvent A was water with 0.04% acetic acid and solvent B was acetonitrile with 0.04% acetic acid. The mobile phase was consisted of solvents A and B. A gradient program was executed in the following steps at 0, 11, 12.1 min: 95:5 (v/v), 5:95 (v/v), 95:5 (v/v) at 12.1 min. The flow rate was 0.4 mL/min, the temperature was 40 °C, and the injection volume was 2 μL. The effluent was alternatively connected to an ESI Q-TRAP (AB Sciex QTRAP 6500 System, AB SCIEX Pet. Ltd, Framingham, MA, USA) (Zhang, Lu & Liu, 2022). The parameters of ESI source operation were set as follows: ion spray voltage (IS, 5,500 V), ion source (turbo spray), source temperature (500 °C), ion source gas I (GSI, 55 psi), gas II (GSII, 60 psi), and curtain gas (CUR, 25 psi), and the collision gas (high). Instrument tuning and mass calibration were performed with 10 μmol/L polypropylene glycol solutions in QQQ and 100 μmol/L in LIT modes. Multiple reaction monitoring (MRM) experiments with collision gas (nitrogen) at 5 psi were performed to obtain QQQ scans (Li et al., 2017). The declustering potential (DP) and collision energy (CE) for individual MRM transitions were measured, then DP and CE were further optimized. For each period, MRM transitions were monitored on the metabolites in elution (Dong et al., 2019). Metabolites were annotated according to the information from the public metabolite database and the MetWare database (MetWare Biotechnology Co., Ltd., Wuhan, China). All mass spectrum peaks were subjected to area integration, which represents the relative content of the corresponding metabolites (Zhang, Lu & Liu, 2022).

Fisher least significant difference (LSD) post hoc test was used to compare the amino acids content at different developmental stages. Variable importance in projection (VIP) ≥ 1 and log₂ |fold change| ≥ 1 were set as the criteria to screen the significant differentially accumulated amino acids. Each sample was performed in triplicate. Principal component analysis (PCA) of free amino acids and quality control samples was conducted by Origin (version 2018).

Transcriptomic analysis

Data were collected as previously described in Lin et al. (2020). Specifically, total RNA was extracted from frozen leaf samples using the Total RNA Extractor (TRIzol) kit (B511311; Sangon, Shanghai, China). The quality and concentration of RNA were determined using a NanoDrop 2000 (Thermo Fisher, Waltham, MA, USA). Two micrograms of the RNA was used for library construction by the VAHTSTM mRNA-seq v2 Library Prep Kit (Illumina, San Diego, CA, USA). First-strand cDNA was synthesized by M-MuLV reverse transcriptase using random hexamer primers, while second-strand cDNA was synthesized by DNA polymerase I and RNase H. The left overhangs were turned into blunt ends by exonuclease/polymerase. AMPure XP system was used to purify library fragments (Beckman Coulter Company, Beverly, MA, USA). PCRs were performed with Phusion High-Fidelity DNA polymerase (Thermo Fisher Scientific Inc., Waltham, MA, USA), universal PCR primers, and Index (X) Primer. PCR products were purified (AMPure XP system), and the quality of the library was determined using a Bioanalyzer 2100 system (Agilent Technologies Inc., Santa Clara, CA, USA). Paired-end sequencing of these libraries was performed on HiSeq X Ten sequencers (Illumina, San Diego, CA, USA) by Novagen Co., Ltd. (Beijing, China). Raw reads of high-throughput sequencing were obtained in the FASTQ file format. De novo assembly of the clean reads were conducted by Trinity (version 2.0.6) (parameter: min_kmer_cov 2) (Trinity Technologies, Irvine, CA, USA). To decrease redundancy, we clustered the transcripts with a minimum length of 200 bp. The longest sequence was preserved and designated a unigene for each cluster (Lin et al., 2020). Read alignment statistics and sample quality features were calculated with SAMtools and RSeQC (version 2.6.1) (Lin et al., 2020).

BLAST program was used to annotate the functions of unigenes against protein databases, including NCBI Nr (NCBI non-redundant protein sequences), Swiss-Prot protein, TrEMBL, ConserGene Ontology (GO), CDD (Conserved Domain Database), Pfam, and KOG (eukaryotic Orthologous Groups) database (e-value < 1 ×1 0⁻⁵). GO (Gene Ontology Database) function annotation information was obtained according to the Uniprot annotation results based on the annotation results of Swiss-+Prot and TrEMBL database. KEGG annotation was conducted by KAAS (KEGG automatic annotation server) (version 2.1). All the assembled unigenes were classified into KEGG pathways using BLASTX against the KEGG database (Zhang et al., 2018). Transcripts per million (TPM) was used to eliminate the influence of gene lengths and sequencing discrepancies (Jiang et al., 2019). DESeq2 (version 1.12.4) was conducted to determine the differentially expressed genes (DEGs) between the two samples at a q-value < 0.001 and log₂ |fold change| ≥ 1 (Lin et al., 2020). Each sample was performed in triplicate. PCA of DEGs was conducted by Origin (version 2018).

Reverse transcription-quantitative PCR (RT-qPCR)

Sixteen candidate genes involved in glycolysis, the tricarboxylic acid cycle (TCA) and amino acid anabolism were chosen and RT-qPCR was conducted to verify the gene expression. The specific primers were designed using Primer Premier 5.0 (Table S2). Total RNA was reverse-transcribed using HiScript II reverse transcriptase according to the manufacturer’s instructions (Vazyme, Nanjing, China). Real-time PCR experiments were performed using a CFX Connect Real-Time PCR system (Bio-Rad Laboratories Inc., Hercules, CA, USA). β-Actin was used as the internal standard gene (Zhao et al., 2018; Lin et al., 2021). The relative expression of genes was calculated by the 2^−∆∆Ct method (Salvatierra et al., 2010). All RT-qPCRs were performed with three biological replicates. All statistical analyses were conducted by SPSS 17.0 software (SPSS Inc. Chicago, IL USA). The data are displayed as the means ± standard deviations (SD).

Integrated analysis of transcriptomic and amino acids data

The amino acids and transcriptomic data were natural logarithm transformed, and the Pearson correlation coefficient (PCC) between amino acids and transcriptomic data was calculated by the cor function in R (version 3.5.1). |PCC| > 0.8 and p < 0.05 were set as the cutoff criteria to assess the association between amino acids and DEGs in C. paliurus leaves over three different developmental stages. Finally, the network was visualized through CytoScape (version 3.6.1) (Cho et al., 2016).

Prediction of key transcription factors and phylogenetic analysis

The key transcription factors (TFs) involved in amino acid biosynthesis were predicted by PlantTFDB (http://planttfdb.gao-lab.org/) (Jin et al., 2017). The expression data were natural logarithm transformed, and the correlation between differentially expressed TFs and other DEGs was calculated by the Corrplot package in R (version 3.5.1) with PCC > 0.8, and p-value < 0.05 were considered to be significant. A visual network graph was generated with Cytoscape software (version 3.6.1).

The protein sequences of C. paliurus TFs were translated from coding sequences in the transcriptome. The bZIP family protein sequences of Arabidopsis thaliana and rice (Oryza sativa japonica) were downloaded from NCBI. Phylogenetic and molecular evolutionary analyses were conducted using MEGA version 5.05 with pairwise distance and the maximum-likelihood algorithm. The reliability of each tree was established by performing 1,000 bootstrap sampling steps.

Free amino acid profiling in leaves of C. paliurus

Metabolite analysis annotated 694 metabolites in S1, S2 and S3 stages, each of which was analyzed using three biological replicates (Table S3). Of the 694 identified metabolites, organic acids (69, 9.94%), nucleotides and their derivatives (59, 8.50%), flavones (57, 8.21%), amino acid derivatives (53, 7.64%) and flavonols (44, 6.34%) accounted for a large proportion (Fig. S1).

The differences in the relative abundance of free amino acid in the leaves of C. paliurus at different developmental stages are shown in Table 1 and Fig. 2 and Table S4. A total of 18 amino acids were detected. Regarding total amino acids, the relative abundance of amino acids in the leaves at the S1 stage and S2 stage was significantly higher than that at the S3 stage (Table 1).When the threshold for a significant difference was set at variable importance in projection (VIP) ≥ 1.0 and log₂ |fold change| ≥ 1, from the S1 stage to the S2 stage, valine and tryptophan in the leaves of C. paliurus at the S1 stage were 3.6-fold and 43-fold lower, while alanine was 3.4-fold higher, than at the S2 stage; from the S2 stage to the S3 stage, the contents of eleven amino acids decreased significantly (Table 1). PCA results showed that there were significant differences in the amino acid content among the leaves at different developmental stages (Fig. S2).

DEGs involved in the amino acid biosynthesis pathway

Transcriptome sequencing of C. paliurus in the S1, S2 and S3 stages generated average numbers of 52,453,341, 48,977,149, and 46,906,739 raw reads, respectively. In all the libraries, 93.59% of the sequences had quality scores greater than Q30. After filtering out the low-quality reads, 49,806,827, 46,320,477, and 42,536,039 clean reads were obtained from the S1, S2 and S3 stages with GC contents of 49.28%, 51.91%, and 52.63%, respectively. A total of 296,593 unigenes were obtained with an average length of 430.11 bp and N50 of 478 bp (Table S5).

A total of 30 DEGs were involved in the amino acid biosynthesis pathway (KEGG map01230, https://www.genome.jp/entry/map01230, Fig. S3 and Fig. 3, Table S4). PCA results showed that there were significant differences in the DEGs among the leaves at different developmental stages (Fig. S4).

Among ten DEGs related to the glycolysis pathway, pgm, pfkA, gapdh, and eno, which encode phosphoglucomutase, 6-phosphofructokinase 1, glyceraldehyde 3-phosphate dehydrogenase and enolase, respectively, were downregulated, while fba and pgk, which encode fructose-bisphosphate aldolase and phosphoglycerate kinase, respectively, were upregulated.

Among six DEGs related to the TCA cycle, the pgh, cs, icd and scs genes, encoding pyruvate dehydrogenase, citrate synthase, isocitrate dehydrogenase, and succinyl-CoA synthetase, respectively, were downregulated. Among five DEGs related to the pentose phosphate pathway, the rpiA gene, encoding ribose 5-phosphate isomerase A, was upregulated; however, the other three genes, pgd, tkt, and tal, encoding 6-phosphogluconate dehydrogenase, transketolase, and transaldolase, respectively, were downregulated.

Among nine DEGs relate to the amino acids metabolism pathway, the glyA, ltaE, gs, aroC, and tat gene, encoding glycine hydroxymethyltransferase, threonine aldolase, glutamine synthetase, chorismate synthase, and tyrosine aminotransferase, were upregulated; however, the other four gene, hsd, alt, ast, trpE, encoding homoserine dehydrogenase 1, alanine transaminase, aspartate aminotransferase, and anthranilate synthase component I, were downregulated.

Validation of the gene expression profiles by RT-qPCR

To validate the transcription profile obtained by RNA-seq data, 16 genes from starch and amino acid metabolism were selected to design specific primers (Table S2) for RT-qPCR analysis. The results indicated that most of the selected genes (except for the PKAA, ENO and LTAE genes) had similar expression patterns to those identified through RNA-seq analysis (Fig. S5). The overall correlation coefficient (R² = 0.64) between RNA-seq and RT-qPCR data was obtained by linear regression analysis (p < 0.0001) and exhibited a good correlation between these two datasets (Fig. S6, Table S8), indicating that the nine transcriptomic datasets were credible.

Integrated analysis between DEGs and amino acids

The integrated analysis of the log₂-transformed amino acids and transcriptomic data showed that 34 DEGs correlated significantly with 18 amino acids, resulting a total of 110 related pairs (34 positive and 76 negative) between the DEGs and amino acids (Fig. 4, Table S6). Seventeen genes were determined to be involved in involved in the amino acid biosynthesis pathway. Among them, six genes (pgm, pfkA, fba, gapdh, pgk, eno) are related with the glycolysis pathway, one gene (scs) is related with TCA cycle, four genes (pgd, rpiA, tkt, tal) are related with pentose phosphate pathway and six genes (tat, ltaE, hsd, glyA, gs), and ast are related with amino acids metabolism (Fig. 4).

Screening of transcription factors related to amino acid biosynthesis

To clearly determine the unknown putative transcription factors controlling the amino acid biosynthesis of C. paliurus, coexpression analysis was employed between amino acids and differentially expressed transcription factor bZIP genes. A total of 60 differentially expressed bZIP genes were identified via comparison with PlantTFDB (Table S9). A total of 78 related pairs were identified. The visualized network showed that a total of 35 nodes were connected by 78 edges (Fig. 5). Only TRINITY_DN82089_c0_g2 had a positive correlation with 11 amino acids; however, the other bZIP genes were negative for amino acids.

Phylogenetic analysis of bZIP family members in C. paliurus

To further investigate the evolution of the C. paliurus bZIPs, we constructed phylogenetic trees either from the full-length C. paliurus bZIPs or from the bZIP domains alone. As shown in Fig. 6, TRINITY_DN89043_c4_g4 clustered with OsbZIP18 and subsequently clustered with TRINITY_DN82089_c0_g3. TRINITY_DN96834_c0_g2 clustered with AtbZIP52. Both of these bZIPs clustered together and then clustered with TRINITY_DN82089_c0_g2.