Transcriptome analysis of the differences between two kinds of Cassia nomame germplasm resources

Plant materials

C. nomame (T4-1 and T4-2) grows in Xinglong and Kuanchen respectively, Chengde City, Hebei Province, China. C. nomame (T4-1, T4-2) leaves of flowering stage with three replicates were collected and gathered, promptly frozen in liquid nitrogen bottles, and cryopreserved to extract RNA and metabolites. Subsequently, the leaves were placed in an ultra-low-temperature refrigerator.

RNA extraction, library construction, and sequencing

Total RNA was extracted from the samples using TRIzol (Invitrogen, Waltham, MA, USA) and purified utilizing RNeasy column (Qiagen, Hilden, Germany). PolyA-tailed mRNA was isolated from total RNA using oligo(dT) magnetic beads. The purified RNA was then fragmented into approximately 300 bp pieces via ion-mediated shearing. First-strand cDNA was generated from the RNA fragments with random hexamer primers and reverse transcriptase. Subsequently, the second-strand cDNA was produced using the first strand as the template. After library construction, the library fragments were enriched by PCR amplification, and the library was selected according to the fragment size. The library size was 450 bp. The quality of the library was verified by Agilent 2100 Bioanalyzer, and the total and effective concentrations of the library were determined. Based on the library’s effective concentration and the required sequencing depth, pools were created by combining uniquely indexed libraries in specific ratios. This pooled library was diluted uniformly to a 2 nM concentration and then denatured with alkali to generate single-stranded DNA. The resulting libraries were subjected to paired-end sequencing on an Illumina platform.

Transcriptome data filtering and assembly

Data filtering using Cutadapt was employed to remove the 3′-end joint sequence, and reads with an average mass fraction lower than Q20 were eliminated. For transcriptome sequencing without reference genome, Trinity software was used to splice clean reads to obtain transcripts for subsequent analysis (Haas et al., 2013). Trinity is a de novo assembly software for transcriptome assembly, which is based on the principle of De Bruijn Graph (DBG) assembly to splice high-quality sequences. The software consists of three independent software modules with the following workflow: Inchworm: a short-sequence library of K-mer length is constructed using high-quality sequences, and the short sequence is extended by overlap of K-mer-1 length between short sequences to obtain a preliminary spliced contig sequence; Chrysalis: the contig sequences are clustered based on overlap between them, and a Bruijn graph is constructed for each class; and Butterfly: deals with these Bruijn diagrams, and finds the path according to the reads and paired reads in the diagram to obtain the transcript. After completion of stitching, a transcript sequence file in FASTA format can be obtained. The longest transcript under each gene is extracted as the representative sequence of the gene, called unigene.

Transcriptome functional annotation and differential gene expression analysis

The databases used for gene function annotation included NR (NCBI non-redundant protein sequences), Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genome (KEGG), Evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG), Pfam, and SwissProt. By comparing and annotating with the NR library, the similarity between the gene sequences of the target species and related species as well as the functional information of the gene of the target species can be obtained. GO annotation was completed using BLAST2GO software, and the default parameters of BLAST2GO were used for annotation (Conesa & Gotz, 2008). The GO annotation results were mapped to GO terms, and the number of genes annotated to the second-level classification was counted. KO and Pathway annotations were mainly performed using the KOBAS annotation system. After completing the KO annotation, the KO was mapped to the corresponding KEGG pathway. The eggNOG comparison annotation was performed on the gene, and the eggNOG number of the best comparison result was assigned to the corresponding gene (Cantalapiedra et al., 2021). Furthermore, each gene was classified into the eggNOG classification directory by using the correspondence between the eggNOG number and eggNOG classification directory.

Transcriptome expression quantification was performed using the RSEM software. With the transcript sequence as a reference, the clean reads of each sample were compared to the reference sequence (Li & Dewey, 2011). Then, the number of reads on each gene was calculated for each sample, and the FPKM value for each gene was computed. We used DESeq2 to analyze the difference in gene expression (Love, Huber & Anders, 2014), and the conditions for screening differentially expressed genes (DEGs) were as follows: | log2FoldChange | > 1, FDR < 0.05. In total, six samples were included and are evenly separated into two groups, with three biological replicates on each group. The statistical power of this experimental design, calculated in RNASeqPower is 0.85 (FDR < 0.05) (Hart et al., 2013).

Functional and pathway enrichment analysis

GO term enrichment analysis was conducted with the topGO package. Statistically significant terms were identified by applying a hypergeometric test to the differentially expressed genes (DEGs) annotated under each term, using the whole genome as background. A false discovery rate (FDR) threshold of <0.05 was applied to define significantly enriched terms. The results were categorized into biological process (BP), molecular function (MF), and cellular component (CC). For each category, the ten most significant terms (lowest P-value) were selected for visualization. The KEGG database systematically analyzes gene functions and links genomic information and functional information (Kanehisa et al., 2023). The KEGG metabolic pathway enrichment analysis was performed using KOBAS online software, and KEGG pathway with FDR < 0.05 was selected as the enrichment result (Bu et al., 2021).

Transcription factor analysis and real-time quantitative PCR

The process of eukaryotic transcription initiation is very complex and often requires the assistance of a variety of protein factors. The transcription factors (TFs) are a class of protein molecules that can specifically bind to specific sequences upstream of the 5’ end of the gene and form a transcription initiation complex with RNA polymerase II to participate in the process of transcription initiation. The transcription factor (TF) prediction was compared with PlantTFDB (Plant Transcription Factor Database) database to predict the TFs and obtain the family information of the TFs (Jin et al., 2017). The DEGs predicted as TFs were statistically analyzed, and according to the family information of the TFs, the number of differentially expressed TFs in each TF family in the comparison group was statistically analyzed.

Using the Total RNA Purification Kit (Promega, Madison, WI, USA), the samples’ total RNA was extracted. DNase I was then used to degrade the RNA and the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA) was utilized to synthesize cDNA. Three biological replicates of qRT-PCR were run using the SYBR green PCR Master Mix. For equivalent loading, the 18S RNA was utilized as an internal standard. The 2^−ΔΔCt analysis method was used to determine the relative expression levels. The SD of three biological replicates was calculated using GraphPad Prism 8 software, and the student’s t-test analysis was performed. List of primers utilized in this investigation was listed in Table S1.

Protein-protein interaction network analysis

Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) is a protein–protein interaction (PPI) database developed by EMBL, and includes the most powerful experimental evidence, data mining, and homologous prediction of PPIs (Szklarczyk et al., 2023). By including the PPI information of the species in the STRING database and based on the results of DEGs analysis, we screened PPI pairs with DEGs as both ends of the nodes and score > 0.95 in the direct database. When there was no PPI information for the species in the STRING database, we selected similar species to compare the sequences and then obtained the relationship between the proteins of the species.

Yeast two-hybrid assay

The full-length coding sequences of DN82231_c0_g1 and DN55_c6_g1 were amplified by PCR and combined with pGADT7 and pGBKT7 vectors, respectively. GAL4 two-line hybrid system was used for Y2H detection assay and the recombinant plasmid was transformed into Y2H golden yeast strain. The transformed yeast cells were cultured at 30 °C for 3 days in selective medium lacking leucine, tryptophan and selective medium lacking leucine, tryptophan, histidine and adenine, respectively. Then, the growth of yeast cells was observed and reflected the interaction between the two proteins.

Transcriptome data assembly and annotation

To explore the differences among various C. nomame germplasm resources, we selected two C. nomame varieties with significant differences in growth and development, stress response and secondary metabolites based on previous experiments (Fig. S1). Then, we performed de novo transcriptome sequencing for the leaves of two C. nomame cultivars which have different growth status and anthraquinone content (Figs. 1A and S1). Based on the transcriptome data after quality control and screening, the clean reads of the six samples were in the range of 5,747,748,258–6,985,352,110 bp, the percentage of Q30 bases was 93.09–93.83%, and the GC content was 40.17% (Tables 1 and S2). The overall sequencing data filtering quality was good, and could be used for subsequent transcriptome analysis. Through transcriptome assembly, we obtained 56,136 independent genes, with length ranging from 302 to 16,566 bp, average length of 1,088 bp (Fig. 1B), and total length of 66,667,714 bp (Table 1). According to the distribution of different lengths of independent genes after assembly, the gene length of 0–500 bp accounted for about 60% of the total length. Furthermore, the length of N50 was 2,096 bp, which indicated that the assembly results were reliable.

To explore the functions of unigenes, we used five databases (NR, KEGG, GO, eggNOG, SwissProt) to annotate the obtained unigenes (Fig. 1C). A total of 34,783 unigenes with annotation information were obtained, 21,353 unigenes were not annotated, and the annotation ratio was 62%. Among them, 32,444, 19,625, 12,380, 30,361, 20,215, and 26,647 unigenes were annotated in the NR, GO, KEGG, eggNOG, Pfam, and SwissProt databases, respectively, with 7,490 unigenes annotated in all the five databases.

Analysis of DEGs in transcriptome

To ensure accuracy of subsequent analysis, we first corrected the sequencing depth and then adjusted the length of the gene or transcript to obtain the RPKM value of the gene. The expression distribution was calculated according to the RPKM value of each gene, which revealed relatively uniform gene expression among different samples, indicating that these results can be used for subsequent DEGs analysis (Fig. 2A). Principal component analysis (PCA) and Pearson correlation coefficient analysis were employed to cluster the samples according to the expression level. The PCA results showed that the three biological replicates of different C. nomame cultivars were distributed in the same region, indicating data integrity and reliability, as well as diverse expression levels of different genes in the two cultivars. Pearson correlation coefficient analysis revealed that the correlation between the two samples was >0.99 (Fig. S2), implying that the samples had good repeatability. Based on the gene expression level in the transcriptome data, we selected genes with differential expression fold | log2FoldChange | > 1 and P < 0.05 as DEGs. Finally, a total of 4,696 DEGs were identified between the two C. nomame data, of which 2,391 and 2,305 DEGs were upregulated and downregulated, respectively (Figs. 2C and S3).

Functional analysis of DEGs

GO is an internationally standardized classification of gene functions, which provides a biological basis for the global characterization of de novo assembled transcripts (Gene Ontology Consortium, 2015). To elucidate the biological function of DEGs in C. nomame, we performed GO functional enrichment analysis of selected 4,696 DEGs (FDR ≤ 0.05), which included three aspects of biology, namely, BP, MF, and CC. The results revealed that the DEGs of the two C. nomame cultivars were mainly enriched in BP. In the BP classification, the GO terms were mainly enriched in the regulation of DNA binding, secondary metabolites, photosynthesis redox, cell division response, pigment biosynthesis, polysaccharide metabolism, hemicellulose metabolism, response to plant hormones, tetrapyrrole biosynthesis process, flavonoid and terpenoid anabolism, gibberellin-mediated signaling pathway, phylloquinone biosynthetic process, quinone biosynthetic process, and other biological processes (Fig. 3 and Table S3). In the MF classification, the GO terms were significantly enriched in TF activity, tetrapyrrole binding, pigment binding, quinone binding, terpenoid synthase activity, gibberellin oxidase activity, sequence-specific DNA binding, rRNA binding, and dioxygenase activity (Fig. 4 and Table S3). In the CC classification, the GO terms were mainly enriched in chloroplast thylakoid, thylakoid photosystem, cell wall plastid, ribosome, photosystem I reaction center, plant cell wall, ATPase complex, etc. (Fig. 5 and Table S3). These results suggested that the differences in biological functions may lead to variations in the growth status and secondary metabolites of different C. nomame cultivars.

KEGG annotation provides information of transcripts related to metabolic process and functions in cellular processes (Kanehisa et al., 2023). We used the KEGG metabolic pathway enrichment analysis to determine the key metabolic pathways of the DEGs of the two C. nomame cultivars. The results revealed that the DEGs of the two C. nomame cultivars were significantly enriched in ribosome, photosynthesis, flavonoid biosynthesis, carotenoid biosynthesis, phenylpropanoid biosynthesis, anthocyanin biosynthesis, nitrogen metabolism, circadian rhythm, α-linolenic acid metabolism, tryptophan metabolism, terpenoid biosynthesis, plant hormone signal transduction, and other metabolic pathways (Fig. 6 and Table S3). KEGG functional enrichment analysis showed that the DEGs played an important role in the biosynthesis and metabolism of secondary metabolites, transcriptional regulation, and other processes.

Transcription factor identification and differential expression analysis

In eukaryotes, TFs are proteins that bind to sequence-specific DNA. They control gene expression by binding to specific DNA sequences to ensure that they are expressed in the appropriate number of cells at the right time throughout the life cycle of cells and organisms (Lambert et al., 2018). TFs are mainly involved in the initial stage of RNA transcription and are the key factors in regulating the gene expression levels (Strader, Weijers & Wagner, 2022). GO enrichment analysis showed that transcription regulator activity, regulation of DNA binding, and other GO terms related to TFs were significantly enriched, indicating that the regulation of TFs in the two C. nomame cultivars played an important role in the formation of their phenotypes. Furthermore, we predicted the TFs of all unigenes and identified 7,309 candidate TFs. The global TF classification noted that the identified candidate TFs belonged to 57 TF families, such as bHLH, MYB, NAC, ERF, C2H2, MYB-related, WRKY, bZIP, C3H, etc. (Fig. 7A).

In addition, we also performed differential expression analysis of TF families, and found that 48 TF families were upregulated or downregulated in the two C. nomame cultivars. Some key regulatory gene families involved in response to abiotic and biotic stresses, such as WRKY, bHLH, MYB, NAC, bZIP, C2H2, HSF, and C2C2-YABBY, were upregulated or downregulated in both the cultivars (Fig. 7B). However, several TF family genes, such as ERF, C2H2, NF-YA, and TALE, were upregulated in T4-2, while TF families such as GRAS, C3H, BBR-BPC, and STAT were highly upregulated in T4-1, indicating that these TFs may be highly significant for the generation of various characteristics of C. nomame. To verify the authenticity of the transcriptome data, we used qRT-PCR technology to confirm the expression levels of individual genes in some TF families (Fig. 8). The results showed that the expression trend of these genes was consistent with RNA-Seq, which further validated our speculation.

Identification and verification of PPI network

Proteins are the main catalysts, structural elements, signal transmitters, and molecular machines of biological tissues. PPI is very important in the coordination of intracellular events, and is the basis of multiple signal transduction pathways in cells and various transcriptional regulatory networks (Von Mering et al., 2002). To describe the relationship between proteins encoded by the DEGs of the two C. nomame cultivars, we used STRING software to construct a PPI network to evaluate the relationship between protein-encoded DEGs. Based on the results of differential gene expression analysis, we screened PPI pairs with DEGs at both the ends and score > 0.95 in the direct database, and identified 675 interactions between proteins encoded by 171 DEGs (Fig. 9A). Most of these interacting proteins were noted to be encoded by upregulated genes, indicating that these proteins may play an important role in the growth and development of C. nomame. Subsequently, we randomly selected a pair of interacting proteins DN82231_c0_g1 (Malectin domain kinesin) and DN55_c6_g1(Gamma subunit of Arabidopsis chloroplast ATP synthase), and verified them by yeast two-hybrid technology. The findings revealed an interaction between the two proteins (Fig. 9B), which further confirmed the authenticity of the interaction network.