Molecular analysis of the 14-3-3 genes in Panax ginseng and their responses to heat stress

Gene family identification and physicochemical property analysis

Genome data for ginseng were generated from the leaves of a four-year-old Korean cultivar ‘Chunpoong,’ and three individuals were used to reduce heterogeneity. The reads were obtained by the Illumina platform (HiSeq2000 and MiSeq), and the sequences were assembled based on paired-end and mate-paired raw data. Genome and annotation files were available from the Ginseng Genome Database (http://ginsengdb.snu.ac.kr/; Jayakodi et al., 2018). The HMM of 14-3-3 genes was retrieved and downloaded from the Pfam database (Bateman et al., 2004) with access number PF00244, and the E-value threshold was 10⁻⁵. The conserved structural domains of the genes were further identified by the SMART (Letunic, Doerks & Bork, 2015) database (http://smart.embl.de/) and NCBI-Conserved Domain Database (NCBI-CDD), and genes lacking the target structural domains were excluded. The ExPASy program (http://www.ExPASy.org/tools; Rueda et al., 2015) was used to predict the features of 14-3-3 proteins, including molecular weight, isoelectric point, stability, and hydrophilicity. Potential signal peptide cleavage sites in proteins were predicted by SiginalP-5.0 (https://services.healthtech.dtu.dk/service.php?SignalP-5.0). The prediction of protein transmembrane helix structure was implemented by TMHMM-2.0 (https://services.healthtech.dtu.dk/service.php?TMHMM-2.0). SOPMA (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa%20_sopma.html) and SWISS-MODEL (https://swissmodel.expasy.org) were used to predict protein secondary and tertiary structures, respectively (Biasini et al., 2014). Finally, the subcellular localization prediction tools WoLF PSORT (https://wolfpsort.hgc.jp) and CELLO (http://cello.life.nctu.edu.tw/) were used to preliminarily determine the location of the family protein or its gene expression product in the cell (Yu, Lin & Hwang, 2004; Horton et al., 2007).

Gene structure and motif analysis

Based on the genome annotation information, the Gff3 file of the 14-3-3 gene family was obtained to analyze the gene structure. The motif was analyzed using the MEME online tool (Bailey et al., 2009; https://meme-suite.org/meme/). The number of motifs was set to ten and the motif width option to 6–50, and the remaining parameters were left on default. Based on TBtools software (Chen et al., 2020), the gene structure and MEME output information were visualized.

Cis-acting elements and GO analysis

The 2,000-bp sequence upstream of the gene coding sequence (CDS) was extracted from the ginseng genomic data by annotation files, and the upstream sequence of the ginseng 14-3-3 genes was extracted from the collection based on the gene ID and uploaded to the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to predict the type and function of the cis-acting elements (Lescot et al., 2002). To further explore the function of 14-3-3 genes, the protein sequences were uploaded to the eggNOG database (http://eggnog-mapper.embl.de/) to obtain GO terms and corresponding descriptive information of the target genes (Cantalapiedra et al., 2021). Each GO term is at a node; the distance from that node to the root node is the level number. Based on the TBtools (GO Level Count) and GO background files (http://geneontology.org/docs/download-ontology/#go_basic), the predicted results were filtered and analyzed at level 2. Results for GO analysis were then visualized according to the number of genes involved in the biological process (BP), cellular component (CC), and molecular function (MF) in the 14-3-3 gene family.

Evolutionary tree analysis of the ginseng 14-3-3 gene family

Daucus carota, A. thaliana, Gossypium hirsutum, and O. sativa 14-3-3 genes were selected to construct an evolutionary tree with the P. ginseng 14-3-3 gene family. D. carrot 14-3-3 genes were retrieved from the UniProt (Apweiler, 2008) protein database (https://www.uniprot.org). Following gene ID, the A. thaliana 14-3-3 genes were obtained from the Arabidopsis Information Resource (https://www.arabidopsis.org/). The G. hirsutum 14-3-3 genes were downloaded from the Cottongen database (https://www.cottongen.org/), and the O. sativa 14-3-3 genes were obtained from the NCBI by accession number. Multiple 14-3-3 protein sequences were aligned by Multiple Sequence Comparison by Log-Expectation (MUSCLE) within 1,000 bootstrap replications (Manuel, 2013). The comparison results were further cut and processed by trimAl, and the evolutionary tree was constructed by IQ-tree (v1.6.12) to filter the most suitable amino acid replacement model. The output tree file was uploaded to the ITOL tool (https://itol.embl.de/itol.cgi) to annotate and optimize the tree diagram (Letunic & Bork, 2021).

Interaction protein and transcription factor regulatory network analysis

The String database (https://cn.string-db.org/) was used to screen for proteins with potential interactions with ginseng 14-3-3 genes (Szklarczyk et al., 2021). Protein sequences were mapped to the most homologous genes, the interaction score was 0.7, and the max number beyond the level 1 shell was no more than 20 interactors. Based on the genome and annotation information, the 600 bp sequence upstream of 14-3-3 genes was extracted by TBtools (Chen et al., 2020). The regulation prediction module of the Plant Transcriptional Regulatory Map database (http://plantregmap.gao-lab.org/) was used to predict transcription factors that have potential regulatory interactions with the input sequences (Tian et al., 2020). The analysis results were downloaded locally, then the protein interaction network and transcription factor regulatory network were constructed by Cytoscape (Shannon et al., 2003).

Analysis of 14-3-3 family expression among different tissue of ginseng

In order to clarify whether the 14-3-3 gene family has tissue or organ expression differences in ginseng, a transcriptomic data analysis was performed (Wang et al., 2015). Raw data for different tissues of ginseng were downloaded from NCBI with the accession number PRJNA302556. Raw data quality was checked with FastQC, and low-quality reads were removed with Trimmomatic. Reads were matched to the corresponding positions in the ginseng reference genome through Tophat v.2.0.10, and the transcriptome was assembled using Cufflinks (Trapnell et al., 2012). Transcripts were quantified using Features Counts, and a matrix was made to measure gene expressions as normalized transcript per million (TPM) values. Next, the TPM values were transformed and normalized in log₂ (TPM+1), and a heat map of the expression of 14-3-3 genes in different tissues was generated with the euclidean and complete parameters for clustering.

RNA extraction and expression analysis of ginseng 14-3-3 genes under heat stress

The relative expression of target genes was probed by qRT-PCR to determine whether ginseng 14-3-3 genes responded to high temperature. Their expression was analyzed under different processing times. Ginseng tissue-cultured ‘Fu Xing 1’ seedlings in their leaf-expansion stage (15–20 d) were subjected to high-temperature stress at 35 °C for 0, 12, 24, and 36 h, and three biological replicates were used. The petioles and leaves of the seedlings were quickly placed in liquid nitrogen at the corresponding stress points and stored at −80 °C. Mixed RNA samples from petioles and leaves were extracted using RNAprep Pure Plant Kit (DP441; TIANGEN, Beijing, China) and quantified by NanoPhotometer N50 (Implen, Munich, Germany), and the integrity was checked on a 1% agarose gel. The mRNA was reverse transcribed into cDNA using the TransScript^® II One-Step RT-PCR SuperMix (AH411; TransGen, Beijing, China). Specific primers were designed and synthesized for quantitative gene analysis using the β-Actin gene (NCBI: No. AY907207) as a reference; the primer sequences are displayed in Table S1. The gene analysis was done by a two-step method based on the TransStart^® Green qPCR SuperMix (AQ101-01; TransGen, Beijing, China) in LightCycler 96 (Roche, Basel, Switzerland) with a reaction system of pre-denaturation at 95 °C for 30 s, and 40 thermocycles of 94 °C for 5 s and 60 °C for 30 s. The obtained melting curve helped determine the specificity of primers and genes. Two sets of technical replicates were set up for each sample. Experimental data were processed and analyzed using the 2^−∆∆CT method (Livak & Schmittgen, 2001) by Cq value. The expression of the target gene in the control group (0 h) was defined as 1, and thus calculate the fold change of gene expression after heat stress. Statistically significant differences were analyzed and are indicated with asterisks with *p < 0.05 and **p < 0.01 based on independent Student’s T-tests at the top of the columns, with error bars representing the SD for three biologic replicates and two technical replicates.

14-3-3 gene family identification and physicochemical property analysis

According to the HMM model combined with SMART and NCBI-CDD co-identification, 42 14-3-3 genes were identified in ginseng with the conserved structural domain of 14-3-3 proteins and named PgGF14-1 to PgGF14-42 (Fig. S1). Their basic characteristics, physicochemical properties, and subcellular localization were analyzed with the support of bioinformatics (Tables 1 and S2). The full lengths ranged from 402 bp for PgGF14-1 to 15,261 bp for PgGF14-42. The CDS and corresponding amino acids ranged from 291 bp and 96 aa for PgGF14-12 to 2,052 bp and 683 aa for PgGF14-42. However, the CDS length of most genes was below 1,000 bp. The proteins’ molecular weight (MW) ranged from 11.27 kDa for PgGF14-12 to 76.98 kDa for PgGF14-42. Except for PgGF14-38, 41, and 42, which had larger molecular weights, the rest of the PgGF14 proteins ranged from 11.27 to 31.35 kDa. All proteins’ isoelectric points (pI) were lower than 7, except for PgGF14-5 and PgGF14-15, indicating that PgGF14s are mainly acidic. Most proteins had an instability index (II) higher than 40, and their GRAVY values were all negative, indicating that PgGF14s are mostly unstable and hydrophilic. The subcellular localization results obtained from WoLF PSORT and CELLO were consistent, indicating that most PgGF14s are cytoplasmic and nuclear, PgGF14-5 and PgGF14-15 are mitochondrial, and only PgGF14-27 is located in the plasma membrane, suggesting a diversity of action sites and functions (Table S2). In addition, none of the proteins had signal peptide shear sites and transmembrane structures, except PgGF14-3 for the latter (Table S2). The secondary structures (Fig. S2) of PgGF14s indicated that most proteins were α-helices, with less β-folding and irregular coiling. Their tertiary structure (Fig. S3) demonstrated that the proteins primarily existed as dimers, and the monomers were mainly composed of multiple reverse parallel α-helices. The protein conformation mostly displayed a cup-like structure with grooves, which may be crucial in protein phosphorylation and target-protein interactions.

Gene structure and motif analysis

The gene intron-exon structure is unique to eukaryotes. It has been linked to the evolutionary process of gene families. The variability of intron-exon structures among different genes of a family also predicts the diversity of gene expansion and species evolution (Xu et al., 2019). The protein sequences of 42 genes were aligned using the MUSCLE program of MEGA 7.0, and the evolutionary tree of PgGF14s was constructed by the maximum likelihood (ML) method (Fig. 1). The alignment results indicated that PgGF14s displayed a high sequence similarity and strong conservativeness between them, but their C-terminus was variable (Fig. S4). Their structure (Fig. 1) also displayed significant differences, and most genes had two to seven exons. Based on their intron-exon structure, the genes were divided into two major classes, ε and non-ε groups. The ε class genes were located in the same major evolutionary branch. The number of exons of non-ε genes was mainly 2–4, except PgGF14-1 with only one exon, PgGF14-27 with six exons, and PgGF14-38 with five exons. The number of exons of ε genes mainly was 5–7, except PgGF14-41 with 16 exons and PgGF14-42 with 17 exons. Genes displaying a similar structure can have variable intron lengths. The analysis of gene motifs (Fig. 1) indicated that PgGF14s were highly conserved, especially their motifs within the same evolutionary branch with good consistency. In contrast, the number and types of motifs in different evolutionary branches differed considerably. The motif sequences are provided in Fig. S5. Motifs 1, 2, and 3 were the most widely distributed in the gene family. Twenty-five PgGF14s contained all motifs, and PgGF14-41, PgGF14-42 from the ε group had two motif 1. Seven genes included only motifs 1, 2, and 3. Other genes varied widely in the number and type of motifs. PgGF14-5 had only motif 3 and motif 8, and PgGF14-15 had motifs 2, 4, and 8. Except for PgGF14-24 and PgGF14-32, nine genes of class ε had motifs 1 to 10. Although another 16 genes in the non-ε group were highly consistent with the motifs of these nine class ε genes, their structures were very different, possibly related to gene duplication events.

14-3-3 gene family cis-acting elements and GO analysis

Gene expression is the result of multi-level and multi-factor regulation. The cis-acting element is a segment of the DNA sequence that affects gene expression and is essential in site binding. The results and details of the PlantCARE analysis about PgGF14s cis-element types and roles are displayed in Table S3. Twenty-two kinds of cis-acting elements were predicted in the upstream sequences of PgGF14s, mainly including hormone-responsive, light-responsive, stress-responsive, and growth and development regulatory elements (Fig. 2). Five hormone elements responding to ABA, methyl jasmonate (MeJA), auxin (IAA), GA, and salicylic acid (SA) were more widely distributed in the upstream regulatory regions of PgGF14s with corresponding ABRR, TGACG, TGA-motif, P-box, and TCA-element, suggesting that multiple hormones may induce the expression of the ginseng 14-3-3 genes. In addition, 32 genes had the anaerobic induction element ARE, and 19 PgGF14s had TC-rich repeats and binding sites for defense and stress response elements. Seventeen PgGF14s had LTR-based low-temperature response elements upstream of the gene. Upstream sequences of ten genes also contained AT-rich DNA binding protein (ATBP-1), the binding site for AT-rich elements. Some genes that contained upstream binding sites for MYB and MYBHv1 were involved in light responses, drought-inducibility, and flavonoid biosynthetic gene regulation. Furthermore, circadian control, endosperm expression, meristem expression, and zein metabolism regulation elements were also found upstream of the 14-3-3 genes. Two elements were specific to certain genes: a seed-specific regulatory element was found upstream of PgGF14-11, and a wound-responsive element upstream of PgGF14-38. The abovementioned results suggest that the 14-3-3 gene family may be involved in physiological processes, such as growth and development regulation, stress response, substance metabolism, and hormone response in ginseng.

GO annotation and enrichment analysis of genes effectively clarify genes’ biological functions and the processes they are involved in. In order to better analyze the results of GO, it is usually necessary to filter all terms. Here, GO terms for PgGF14s were obtained by eggNOG and further enriched by TBtools at level 2 for analysis. Details of GO annotations and analysis results are displayed in Table S4. PgGF14s were enriched in biological process, cellular component, and molecular function process terms (Fig. 3). All members were enriched in cellular developmental process (GO:0048869), cellular metabolic process (GO:0044237), primary metabolic process (GO:0044238), multicellular organismal process (GO:0032501), and cellular and anatomical entity (GO:0110165). Forty PgGF14s were related to responses to abiotic stimuli (GO:0009628) and cellular response to stimuli (GO:0051716), and 30 PgGF14s were found to respond to biotic stimuli (GO:0009607). These findings reflected the involvement of PgGF14s in biological processes in response to stresses from a gene function level. PgGF14s were also associated with biosynthesis and substance metabolism, such as biosynthetic processes (GO:0009058), primary metabolic processes (GO:0044238), secondary metabolic processes (GO:0019748), and the life activities of cells (GO:0016049, GO:0008219). From the perspective of molecular function, PgGF14s were mainly associated with the catalysis of specific reactions and protein-containing. The enrichment results of the cellular component also predicted a wide distribution of 14-3-3 genes in cells. Overall, the results of the GO analysis indicated that ginseng 14-3-3 genes were involved in various life activities, which is consistent with the functional diversity of the 14-3-3 gene family.

Evolutionary tree analysis of the ginseng 14-3-3 gene family

Fourteen D. carota 14-3-3 genes were identified in the Uniprot database and named DcGF14-1 to 14. The 14-3-3 genes of A. thaliana (AtGRF1–13) (DeLille, Sehnke & Ferl, 2001), G. hirsutum (GhGF14-1–33), O. sativa (OsGF14a–g) (Yao et al., 2007), D. carota (DcGF14-1–14), and P. ginseng (PgGF14-1–42) were used to construct a phylogenetic tree with a total of 110 genes (Fig. 4). The gene numbers and corresponding protein sequences are listed in Table S5. Except for PgGF14-2, OsGF14b, OsGF14d, and OsGF14e,106 genes were divided into four major branches, and PgGF14s were distributed in all branches. The phylogenetic tree indicates that most 14-3-3 genes of the selected species belong to the non-ε group, with 79 genes, including 31 of P. ginseng, 25 of G. hirsutum, nine of D. carota, eight of A. thaliana, and six of O. sativa. Gene homology analysis indicated that 34 of 42 PgGF14s existed in pairs with higher sequence similarity (Fig. S6), possibly because gene replication events produced paralogous genes. In the evolutionary tree, P. ginseng and D. carota 14-3-3 genes had the highest homology, indicating they might be orthologous genes with a close genetic relationship.

Interaction protein and transcription factor regulatory network analysis

Sequence alignment revealed that the highest homology with the target gene in the string database was D. carota. The PgGF14s were mapped to eleven 14-3-3 proteins of D. carota, and some gene pairs displayed more than 90% homology (Table S6). The highest homology correspondence was used to construct an interaction network map of ginseng 14-3-3 proteins (Fig. 5, Table S6). Each node represents a different protein, and the linkage between the nodes indicates the interaction between two proteins. The results indicated that the PgGF14s might interact with functional genes, such as pyruvate kinase (PK), ATP-binding cassette transporter (ABC), glycoside hydrolase family 13, LHC, matK, NPH3, cation transport ATPase (P-type) (TC 3.A.3), and members of transcription factors MYB, heat shock factor (HSF), and Basic/Helix-Loop-Helix (bHLH). Some genes were predicted to interact with proteins containing gene-specific structural domains, such as the Hist_deacetyl, Aminotran_1_2, and BES1_N. The interactions between proteins were diverse. One gene may interact with multiple proteins, consistent with the functional diversity of the 14-3-3 genes and the possibility of many regulatory factors determining the gene expression. Moreover, several 14-3-3 genes were predicted to interact with other genes, such as HSF, MYB, and NPH3. Only PgGF14-20 had potential interactions with all three members of the TC 3.A.3 family. The gene families and functions of some proteins that might interact with PgGF14s in the database remain unclear.

The PgGF14s 600 bp promoter sequences were uploaded to the PlantRegMap database, the reference species was D. carota, and the other parameters were left as default. The analysis identified 355 regulations between 121 TFs and 42 ginseng 14-3-3 genes, and 38 TFs were over-represented with p ≤ 0.05 (Table S7). The potential regulatory relationships between 38 TFs and 42 PgGF14s were constructed as a network diagram, and TFs of the same type were categorized (Fig. 6). Thirteen categories of transcription factors with potential regulatory effects on PgGF14s were predicted, among which MYB, bZIP, and Dof genes accounted for a large proportion, all of which are crucial regulatory transcription factors in plants and are involved in plant abiotic stress response processes. In addition, AP2, MIKC-MADS, and BBR-BPC transcription factors may have potential regulatory roles on several PgGF14s.

Transcriptome data analysis of different ginseng tissues

The species used for transcriptome analysis was a four-year-old Jilin ginseng ‘Damaya’, whose following parts were analyzed: fiber root (SRR2952867), leg root (SRR2952868), main root epiderm (SRR2952869), main root cortex (SRR2952870), rhizome (SRR2952871), arm root (SRR2952872), stem (SRR2952873), leaf peduncle (SRR2952874), leaflet pedicel (SRR2952875), leaf blade (SRR2952876), fruit peduncle (SRR2952877), fruit pedicel (SRR2952878), fruit flesh (SRR2952879), and seed (SRR2952880). The raw TPM values of PgGF14s in different tissue sites are displayed in Table S8. A heat map for gene expression was drawn based on the TPM of PgGF14s in log₂ (TPM+1) (Fig. 7). Nine genes were not detected in this sequencing (PgGF14-2, 4, 5, 6, 12, 13, 15, 17, and 24). PgGF14-3 and PgGF14-16 were expressed just in one tissue and with very low gene abundance, so these genes were not represented in the heat map. Among the remaining 31 genes, 26 (83.9%) were expressed in all tissues, and five (16.1%) were expressed in most tissues. The gene abundance of different members varied greatly. For example, PgGF14-25, 30, and 34 were highly expressed, PgGF14-20, 33, and 37 were moderately expressed, and PgGF14-22, 23, and 36 were lowly expressed. Overall, the expression levels of PgGF14s in different organs were variable. The 14-3-3 gene abundance was low in seeds but relatively high in roots, stems, leaves, and fruits. Meanwhile, the expression of 14-3-3 genes displayed tissue specificity. PgGF14-30 displayed a high abundance in different tissues, especially in the fruit flesh. Some genes were expressed at relatively high levels in the arm root, fiber root, leg root, and rhizome, such as PgGF14-7, 8, 9, 11, and 20. PgGF14-28, 30, and 36 were more abundant in the fruit flesh, fruit peduncle, and fruit pedicel.

Expression analysis of 14-3-3 genes under heat stress

All 14-3-3 genes were detected in the experiment, and the melting curves indicated specificity to the primers and genes. Raw data for qRT-PCR and specific p-values are displayed in Table S9. The results of the qRT-PCR analysis indicated that, overall, the expression patterns and levels of PgGF14s differed between the three treatment times at 35 °C. Representative examples of PgGF14s with different expression patterns are displayed in Fig. 8, and the relative expression for the remaining genes are displayed in Fig. 9. Several genes—PgGF14-6, 13, 14, 16, 18, 23, 24, 31, 38, 39 and 42—were downregulated at 12 h of treatment; their expression rebounded at 24 h and decreased again at 36 h. The expression levels of PgGF14-7, 8, 9, 10, 19, 20, 21, and 34 were high at 12 and 24 h and lower at 36 h, displaying a trend of rising first and falling later. Some genes responded significantly to heat at a particular time point. PgGF14-19 was significantly upregulated at 12 h of high-temperature stress, with a 3.759 times higher expression than the control. PgGF14-18 and 22 were remarkably upregulated at 24 h of stress, with 4.635 and 3.629 times higher expression than the control, respectively. The expression levels of PgGF14-12, 25, 26, 28, 29, 30, 35, and 36 were not very affected at 12 and 24 h treatment times but were significantly downregulated at 36 h of stress. In addition, PgGF14-4 was significantly downregulated at different treatment times, and its expression tended to decrease with increased stress time. On the other hand, PgGF14-5 was significantly upregulated at different treatment times. The relative expression of PgGF14-4 and PgGF14-5 also changed considerably at different processing times. The same trend and apparent response at several time points suggest that PgGF14-4 and PgGF14-5 may be continuously involved in the response of ginseng to high-temperature stress and may have a more critical regulatory role. Among the 42 14-3-3 genes in ginseng, 38 responded more significantly to the high-temperature stress, indicating that the ginseng 14-3-3 gene family may be involved in response to abiotic stresses, particularly heat stress. The differences in expression patterns and levels of different members also indicated the diversity and variability of PgGF14s functions.