Illumina paired-end sequencing data and De novo assembly
A total of 97.86 Gb sequencing data were obtained, including 676,944,532 raw reads and 668,794,834 clean reads with the base average error rate below 0.03%. After a rigorous quality check and data filtering, 668,794,834 high-quality clean reads were obtained with 98.36% Q20 and 94.89% Q30 bases. The GC percentage in ground parts (leaves, stem) and underground parts (roots) were an average of 44.6% and 44.65%, respectively (Table 1). The total number of unigenes with paired-end reads was 110,779, and the total length of the unigenes was 92,992,355 bp, with an average length of 839.44 bp and the N50 and E90N50 value of 1,540 and 3,117 bp, respectively. In the 110,779 unigenes, 24,442 unigenes (22.1%) were greater than 1 kb in length (Table 2).
Functional annotation
The transcriptome sequences were annotated by six databases (NR, Swiss-Prot, Pfam, COG, GO and KEGG) to obtain similarity sequence and the corresponding annotation information. Gene annotation showed that 63,863 unigenes were successfully annotated in Pfam, Swiss-Prot, NR, COG, KEGG, GO databases. The number and mapping rates of unigenes against the Pfam, Swissprot, GO, COG, KEGG databases were 42.05%(46,581)ˎ 42.03%(46,559)ˎ 40.05%(45,036)ˎ 48.11%(53,295)ˎ 29.74%(32,950), respectively (Figure 1a)(Table 2). 34,092 unigenes had high similarity (greater than 80%) in mapped sequences with Nr database and 36,632 unigenes (63.39%) had significant homology (e-value < 10−30) (Figure 1b and 1c). Species distribution analysis showed that only 20,060 unigenes (34.71%) had high homology with the genes from Quercus suber, followed by Abrus precatorius (7,272, 12.58%), Spatholobus suberectus (4,505, 7.8%), while 9,478 unigenes had high homology with sequences from other organisms (Figure 1d).
Functional classification
The functions of all unigenes were classified by using the Nr annotation and Gene Ontology (GO) classification, and these terms could be grouped into biological process (60,270), cellular component (75,836), and molecular function (61,742) of the three main categories to annotations (Additional file 2) (Figure 2). In the molecular function group, we found unigenes related to “binding” (26,213, 26.21%) and “catalytic activity” (24,087, 53.48%). For the cellular component category, “cell part” (22,578, 50.13%), “membrane part” (13,816, 30.68%), “organelle” (12,882, 28.60%), “protein-containing complex” (9,262, 20.57%), “membrane” (5,272, 11.71%) represented the majority of unique sequences. Among biological process category, unigenes assigned to “cellular process” (21,820, 50.13%), “metabolic process” (19,117, 42.45%) were the most abundant. A high percentage of genes were grouped into the “biological regulation” (5,564, 12.35%), “response to stimulus” (3,534, 7.85%), “cellular component organization or biogenesis” (3,533, 7.84%) categories, and “localization” (3,185, 7.07%).
To further understand the biological functions and interactions of transcripts, the unigenes of assembled sequences were assigned by the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The result showed that a total of 31,215 unigenes were assigned to 148 KEGG pathways using BLASTx, with an e-value < 1e-5, and were assigned to six main categories. In the six main categories, metabolism was the biggest category (14,751), followed by information processing (9,817), cellular processes (2,483), organismal systems (2,228) and environmental information processing (747) (Additional file 3) (Figure 3). “Translation” had the largest number of unigenes (6,455 unigenes) followed by “Carbohydrate metabolism” (4,069 unigenes), “Energy metabolism” (2,764 unigenes), “Amino acid metabolism” (2,457 unigenes), “Folding, sorting and degradation” (2,036 unigenes), and “Transport and catabolism” (1,853 unigenes). The metabolic pathways in our study were: “carbohydrate metabolism” (4,069 unigenes), “energy metabolism” (2,764 unigenes), “amino acid metabolism” (2,457 unigenes), “lipid metabolism” (1,506 unigenes), “metabolism of other amino acids” (1,123 unigenes), “metabolism of cofactors and vitamins” (789 unigenes), “nucleotide metabolism” (716 unigenes), “Biosynthesis of other secondary metabolites” (630 unigenes), “Metabolism of terpenoids and polyketides” (451 unigenes), “Glycan biosynthesis and metabolism” (242 unigenes). KEGG genetic information processing included “Translation” (6,455 unigenes), followed by“folding, sorting and degradation” ( 2,036 unigenes) and “transcription” (1,054 unigenes). In the environmental information processing category, the most abundant subcategories were “Environmental adaptation” (1,818 unigenes).
Analysis of Differentially Expressed Genes
The unigenes from different tissues of I. stachyodes were compared using assembled data as a reference (Figure 4a). Under the criteria of p-adjust < 0.05 and |log2FC| > = 2, a total of 11,648 differentially expressed genes (DEGs) between leaf and stem were identified. Among them, 9,058 genes were up-regulated, and 26,711 genes were down-regulated. In addition, 35,490 (7,589 up-regulated and 27,901 down-regulated), 37,234 (10,250 up-regulated and 26,984 down-regulated), 12,467 (5,143 regulated and 7,324 down-regulated), 14,989 (7,866 up-regulated and 7,123 down-regulated), 11,648 (7,430 up-regulated and 4,218 down-regulated) were identified in the comparison of leaf vs rootⅠ, leaf vs rootⅡ, stem vs rootⅠ, stem vs rootⅡ, rootⅠ vs rootⅡ, respectively (Figure 4b). To obtain a comprehensive understanding of DEGs, gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG)-based functional enrichment were conducted. According to GO assignments, a total of 27,999 up-regulated DEGs (Additional file 4) and 56,166 down-regulated DEGs (Additional file 5) were divided into three main categories: biological process, cellular component, and molecular function. Overall, the up-regulated and down-regulated DEGs in different goups significantly enriched in the same or different GO terms (Additional file 6). The main GO terms of our research were significantly enriched in the down-regulated DEGs. In group leaf vs rootⅠ or leaf vs rootⅡ, most of those GO terms in group leaf vs rootⅠ are related to metabolism pathways, including “secondary metabolite biosynthetic process”, “Phenylpropanoid biosynthetic process”, “lignin biosynthetic process”, “jasmonic acid metabolic process”, “jasmonic acid biosynthetic process”, “gibberellic acid homeostasis”, “ribosomal small subunit assembly”, “ribonucleoprotein complex biogenesis”. Whereas, in group rootⅠ vs rootⅡ, the down-regulated DEGs were mainly enriched in terms related to “triterpenoid metabolic process”, “triterpenoid biosynthetic process”, “plant-type secondary cell wall biogenesis”, “plant-type cell wall organization or biogenesis”, “plant-type cell wall biogenesis”, “Phenylpropanoid metabolic process”, “Phenylpropanoid catabolic process”, “lignin catabolic process, external encapsulating structure organization, cell wall biogenesis”.
Among the KEGG pathway analysis, biosynthesis of secondary metabolites such as “Phenylpropanoid biosynthesis” and “Flavonoid biosynthesis” represented the top twenty enriched KEGG pathways especially in the up-regulated DEGs of group leaf vs stem, leaf vs rootⅠ, stem vs rootⅠ, rootⅠ vs rootⅡ, and the down-regulated DEGs of group stem vs root and rootⅠ vs rootⅡ (Additional file 7, Additional file 8). Notably, the down-regulated DEGs in group leaf vs rootⅠ and leaf vs rootⅡ were significantly enriched in Ribosome, and only group stem vs rootⅠ had up-regulated DEGs that significantly enriched in the “Phenylalanine metabolism pathway”, containing 70 upregulated DEGs (Figure 4c). In conclusion, DEGs analysis shows that the PAs accumulation in I. stachyodes root might be due to differentially expressed genes involved in these biological metabolic processes.
Expression patterns of PA biosynthesis potential pathway structural genes in different organs
A large number of flavonoids and PAs were detected in I. stachyodes roots, the two-year-old I. stachyodes root with red color usually used as the medicinal harvesting standards, which means the content of flavonoids in the root of two-year-old I. stachyodes is higher[27]. To further investigate these important findings, the transcriptome of different tissues and roots with different growth years (leaf, stem, rootⅠ, and rootⅡ) were compared to dig out the key genes in the metabolism of red root related to growth years of I. stachyodes. In total, 75 unigenes that encoded 14 enzymes in the flavonoid and PA biosynthesis pathways were identified using BLASTp with previously identified Arabidopsis thaliana genes annotated in KEGG and additional literature[37-39] (Additional file 9). The normalized expression profiles of all the putative PA biosynthesis unigenes found in the I. stachyodes transcriptome were shown in Figure 5. The biosynthesis pathway structural genes of PAs have been mainly divided into three parts[34]. The general biosynthetic genes (GBGs) including PAL, C4H and 4CL are marked in blue, the EBGs including CHS, CHI, F3H, and F3′H are marked in green, while the LBGs including DFR, ANS, ANR, LAR, UGT72L1, MATE and LAC15 are marked in red.
In general, most of the pathway structural genes had higher expression in rootⅠ tissue (IsPAL7, IsC4H1, IsC4H2, Is4CL11, IsCHS1, IsCHS4, IsCHS7, IsCHS8, IsCHI5, IsUGT72L1-5, IsUGT72L1-10, IsUGT72L1-12) or in rootⅡ tissue (IsF3′H1, IsANR2, IsLAR2, IsUGT72L1-2, IsUGT72L1-3, IsUGT72L1-4, IsUGT72L1-8, IsMATE2, IsMATE3), and had lower expression in rootⅠ tissue (IsCHS6, IsDFR2, IsANS, IsUGT72L1-6) or in rootⅡ tissue (Is4CL1, Is4CL18, IsCHS3, IsCHI4) (Figure 5). Interestingly, both the genes upregulated in rootⅡ and the genes down regulated in rootⅠ belong to EBG and LBG, mostly LBG. We focused more on 9 pathway genes (IsF3′H1, IsANR2, IsLAR2, IsUGT72L1-2, IsUGT72L1-3, IsUGT72L1-4, IsUGT72L1-8, IsMATE2, IsMATE3) which significantly upregulated in rootⅡ, the in-depth study of them can clarify the molecular mechanism of PA accumulation in I. stachyodes roots.
Network analysis
PA biosynthesis is controlled by regulatory networks that consist of TFs or regulatory complexes in different species[13, 40]. In order to comprehensively reveal the regulatory network of PA biosynthesis, the expression data of the PA pathway genes and all the TFs were extracted from I. stachyodes transcriptome dataset (Additional file 9; Additional file 10), and were applied to co-expression analysis using the TGMI algorithm. The triple gene blocks identified by TGMI algorithm with a cut-off significance level of 0.05 (Additional file 11). The interference frequencies of TFs on pathway genes were displayed in the descending order (Additional file 12). Among the top 185 TFs regulators, which interfere with the pathway genes with highest frequencies, in the lists identified by TGMI, 34 IsMYBs, 18 IsbHLHs, 15 IsWRKYs, 9 IsMADSs, 3 IsWIPs are known PA pathway regulators supported by literature. These TFs were further combined to generate a circular network, as shown in the Figure 6, with the TFs arranged in a clockwise direction, from the most frequent to the least frequent. Each directed edge from a TF to a pathway gene represents a regulatory relationship. It is perceivable that the core pathway regulator MYB, highlighted in a light coral color. As shown in the figure 6, MYB does play a central role in the expression regulation of pathway structural genes[13, 14].
In brief, co-expression analysis identified numerous potential interactive regulators of PA biosynthesis, involving 34 IsMYBs, 18 IsbHLHs, 15 IsWRKYs, 9 IsMADSs, and 3 IsWIPs (Additional file 13) were chosen for heatmap analysis (Figure 7). It is worth noting that IsWRKY45, IsMYB24, IsbHLH33, IsMYB80, IsMYB9, IsMYB52, IsMADS12, IsMYB68, IsbHLH14, IsMYB79, IsMYB69, IsMYB23 exhibited a higher expression level in rootⅡ. Genes with same or similar expression patterns are often under the regulation of the same molecular mechanism[41]. Thus, we should focus on the network in rootⅡ tissue in order to further dig out key regulators affecting the synthesis of PAs in roots. Two mainly different subnetworks were detected. In sub-network one, four MYBs (IsMYB23, IsMYB79, IsMYB9, IsMYB80), one IsWRKY45, one IsMADS12, and 2 LBGs (IsMATE2, IsMATE3) were co-expressed. Meanwhile, in sub-network two, four MYBs (IsMYB24, IsMYB52, IsMYB68, IsMYB69), two bHLHs (IsbHLH14, IsbHLH33), and only with 1 LBGs IsANR2 co-expressed.
Identification of key regulators of PA biosynthesis genes in I. stachyodes roots
R2R3-MYB generally plays a central role in regulating target genes in PA pathways [13, 14]. To further screen out the PA-related R2R3-MYB proteins and predict their functions, we constructed a phylogenetic tree comprising the 34 IsMYBs proteins along with 126 Arabidopsis R2R3-MYB proteins and 16 proteins related to this process in other plant species (Figure 8). IsMYB79 with higher expression level in rootⅡ than other tissues, was clustered in subgroup 6 of the MYB gene family, such as AtMYB90, AtMYB75, AtMYB114, AtMYB113 in Arabidopsis thaliana[42]. The overexpression of AtMYB75 or AtMYB90 in purple transgenic tobacco plants strongly enhances anthocyanin contents via upregulating all of the anthocyanin biosynthetic genes[43]. IsMYB24, another higher expressed in rootⅡ, which was clustered in subgroup 5 of the MYB gene family, AtMYB123, and TT2-type genes were involved in anthocyanin and PA biosynthesis regulation [13, 44]. IsMYB75 clustered in subgroup 7 and IsMYB22 clustered in subgroup 5, but they specifically expressed in leaf tissues not in roots (Figure 8).
To determine the characterization of IsMYB24 and IsMYB79, homologous sequence alignment was carried out using deduced amino acid sequences and other published flavonoid-related genes amino acid sequences (Figure 9). The results show that IsMYB79 and IsMYB24 have the general characteristics of R2R3-MYB gene family, and contained R2 and R3 domains (Figure 9). IsMYB79 was closely related to other published anthocyanin-related MYBs, which were promoting pigmentation, such as CmMYB6[45], MaAN2[46], LrMYB15[47], StMYB113[48], PpMYB10[49], MrMYB1[50] had been studied in model plant tobacco, the molecular mechanism of regulating anthocyanin accumulation has been basically clarified. EsMYBA1 influence pigmentation in the leaves, flowers, and flower buds[51]. LhMYB12 and LhSorMYB12 in Lilium species controls anthocyanin pigmentation in whole tepals[52]. The transcriptional activation of RsMYB1[53] resulted the anthocyanin pigmentation. The highly homology indicated that the function of these MYBs were similar. In addition, IsMYB24 was closely related to VvMYBPA2, which play crucial roles in regulating PA biosynthesis[54]. Therefore, IsMYB79 and IsMYB24 were similar to other flavonoid-related genes, which may play an important role in promoting root pigmentation.