Transcriptome sequencing.
Transcriptome sequencing of root samples collected at different time points from two mulberry varieties infected with M. enterolobii was carried out on the Illumina platform (Table 2). The proportion of high-quality sequences after filtration of all samples was > 99%, Q20 (%) was > 98%, and Q30 (%) was > 94%, indicating that the transcriptome data were good quality and the sequencing results highly reliable (Table 4).
Table 1
Unigene ID | Forward primer | Reverse primer |
0014355 | TCTCCCGTCTCATCGAGGCTTA | ACCAAGCACGTTGGCTCTCTA |
0066216 | TCTGGCCAAGCAGATCAACGA | ATCTTTCCCCTGAGCTTGACGC |
0050791 | GTGCTGTTCTTCGACGCTCTCA | TTTCCTCCCGATCCGAACTCCA |
0065349 | AACGTCGCAAAGGGTTCTCC | TGCCACAAGCCCTATTGCAG |
0062337 | TTTCGGTTGGCGCGAACATC | AGCAAGCAACGACGGAACCA |
0023187 | ACCGCCACAGGCTAATGAAGT | TGGTGTCTGTCCCTTGGCTT |
0047645 | TGGCTGATGTTTGTGGCCTTGA | AGTTGCTCTTTGATTGGTGGTGG |
0026072 | TTCTGCTGATACGGCGTCCT | TAGGGCCGCAGAGCTTGATT |
0002892 | GGCATGGCTACTTCCATGTGGT | ATGGTTGGGTTTGCGCTTCC |
0008985 | AGCTCCATTTGGAACGCGGA | GTCGAATCCCTTGTATGAGGCCA |
0059588 | GCAGCAGATGATGGCCTTTC | TGGCGCTTCTTCTCCTTGAG |
0059369 | CGTCAAGGCCAAGACCAATTTCC | CGCCGACGTTGATGTTGTTCTC |
0017931 | ATTGGAATCCGAAGCCGTGGA | TGTCGCCTTTCCGTCCAATCT |
0004573 | GCATCTCTAGTGGCTGCTGCTT | TGGCCGTCTTCAGAGGTTGT |
0018547 | TTTTCCGCCGGGTCGATGTT | TTGACGGCAGTGGCATGAGA |
0026157 | AGAAGGAGGAGGAAGAGGAAGTGT | TGAGGCACGTTTGCGAAGCA |
0070844 | AGTTCTTCATCCGTCGCTCCA | ACAATCTCGCAGACGGCACA |
0026640 | TGATCCAACGAGGGGGTCAT | GGCTTGAATGAAGAAGTTGCTCGG |
0014688 | AGTGGGCGACCAGATCAACA | TGGTTGAGAGTGCCGGTCTT |
0069934 | TCGACATTTTCGGGAAGGCCA | AGCATCATCGGAGTCGCTGT |
0071495 | GCAATGGGGACTGTGACAACA | AGGCATGCCAAAGCTCCAAG |
β-actin | AGCAACTGGGATGACATGGAGA | CGACCACTGGCGTAAAGGGA |
Table 2
Component name | Add amount |
2 × SYBR Green Pro Taq HS Premix | 10µL |
Primer F (10 µM) | 0.4µL |
Primer R (10 µM) | 0.4µL |
Template | 2µL |
RNase free water | up to 20µL |
Table 3
Conditions for qPCR reaction
Reaction stage | Temperature | Time | Cycle |
Predenaturation | 95℃ | 30sec | 1 |
Amplification | 95℃ | 5sec | 40 |
65℃ | 30sec |
Melting curve | 95℃ | 30sec | 1 |
65℃ | 30sec |
95℃ | - |
Condensation | 40℃ | 10sec | - |
Note: 1. When setting the annealing temperature of the amplification stage to 65 °C, you need to select Single in the Acquisition Mode (Acquisition Mode), It means that fluorescence signals will be collected during this period. |
2. It is necessary to select Continuous in the acquisition mode of the second 95 °C stage of the melting curve without setting the time. |
Table 4
Transcriptome sequencing data analysis
Sample | Raw reads, na | Clean reads, n (%)b | Q20 (99% base call accuracy), %c | Q30 (99.9% base call accuracy), %d | GC content, %e |
KB0d-① | 60909044 | 60493992 (99.32%) | 98.2 | 94.71 | 46.43 |
KB0d-② | 62997456 | 62491878 (99.2%) | 98.68 | 96.05 | 45.91 |
KB0d-③ | 66197912 | 65647714 (99.17%) | 98.5 | 95.33 | 46.73 |
GB0d-① | 54258632 | 53912710 (99.36%) | 98.31 | 94.91 | 46.89 |
GB0d-② | 64601656 | 64131504 (99.27%) | 98.78 | 96.17 | 46.9 |
GB0d-③ | 63892702 | 63436502 (99.29%) | 98.83 | 96.35 | 47.0 |
KBjz8d-① | 59244696 | 58881156 (99.39%) | 98.33 | 94.97 | 47.05 |
KBjz8d-② | 60152896 | 59777740 (99.38%) | 98.22 | 94.69 | 47.23 |
KBjz8d-③ | 57259132 | 56906180 (99.38%) | 98.31 | 94.91 | 47.33 |
KBck8d-① | 60755908 | 60391936 (99.4%) | 98.32 | 94.92 | 47.74 |
KBck8d-② | 52040558 | 51686572 (99.32%) | 98.26 | 94.81 | 47.29 |
KBck8d-③ | 64174538 | 63719678 (99.29%) | 98.27 | 94.82 | 47.34 |
GBjz8d-① | 68430676 | 67917050 (99.25%) | 98.22 | 94.72 | 47.96 |
GBjz8d-② | 63719816 | 63292076 (99.33%) | 98.26 | 94.81 | 47.32 |
GBjz8d-③ | 61292316 | 60748664 (99.11%) | 98.22 | 94.77 | 46.92 |
GBck8d-① | 48996092 | 48647074 (99.29%) | 98.28 | 94.88 | 47.62 |
GBck8d-② | 85457146 | 84877082 (99.32%) | 98.29 | 94.9 | 48.08 |
GBck8d-③ | 63423860 | 62959404 (99.27%) | 98.16 | 94.6 | 46.78 |
KBjz17d-① | 81179024 | 80436262 (99.09%) | 98.2 | 94.69 | 47.1 |
KBjz17d-② | 64554492 | 64163368 (99.39%) | 98.33 | 94.93 | 47.82 |
KBjz17d-③ | 50325862 | 50006746 (99.37%) | 98.32 | 94.91 | 47.69 |
KBck17d-① | 57315338 | 56890984 (99.26%) | 98.25 | 94.83 | 48.96 |
KBck17d-② | 58278390 | 57858956 (99.28%) | 98.32 | 95.01 | 49.78 |
KBck17d-③ | 60404646 | 59916916 (99.19%) | 98.85 | 96.37 | 48.53 |
GBjz17d-① | 60396566 | 59987776 (99.32%) | 98.4 | 95.15 | 47.6 |
GBjz17d-② | 73025484 | 72582168 (99.39%) | 98.33 | 94.98 | 47.74 |
GBjz17d-③ | 72193056 | 71733542 (99.36%) | 98.2 | 94.63 | 47.68 |
GBck17d-① | 54270666 | 53902126 (99.32%) | 98.44 | 95.29 | 46.68 |
GBck17d-② | 64310188 | 63735094 (99.11%) | 98.72 | 96.18 | 46.82 |
GBck17d-③ | 55197462 | 54697628 (99.09%) | 98.69 | 96.14 | 47.8 |
KBjz23d-① | 59095200 | 58691120 (99.32%) | 98.31 | 94.93 | 47.14 |
KBjz23d-② | 79210406 | 78762440 (99.43%) | 98.33 | 94.96 | 47.72 |
KBjz23d-③ | 61962352 | 61612726 (99.44%) | 98.41 | 95.16 | 47.66 |
KBck23d-① | 57595244 | 57045606 (99.05%) | 98.78 | 96.32 | 49.15 |
KBck23d-② | 63076062 | 62601422 (99.25%) | 98.77 | 96.3 | 47.09 |
KBck23d-③ | 63769422 | 63342658 (99.33%) | 98.83 | 96.38 | 47.32 |
GBjz23d-① | 64686018 | 64110442 (99.11%) | 98.24 | 94.82 | 47.08 |
GBjz23d-② | 69797772 | 69378514 (99.4%) | 98.32 | 94.94 | 48.46 |
GBjz23d-③ | 47339138 | 47004460 (99.29%) | 98.22 | 94.75 | 49.18 |
GBck23d-① | 52537608 | 52003500 (98.98%) | 98.72 | 96.21 | 49.64 |
GBck23d-② | 60610656 | 60072490 (99.11%) | 98.65 | 95.94 | 47.82 |
GBck23d-③ | 63395176 | 62774902 (99.02%) | 98.71 | 96.13 | 48.32 |
Notes: aNumber of original data reads; bNumber of high-quality sequences after filtering and (percentage of raw reads; cPercentage of bases with an accuracy of ≥ 99% after filtering; dPercentage of bases with an accuracy of ≥ 99.9% after filtering; ePercentage of G/C nucleotides in the sequence. |
Results of de novo assembly.
Trinity software was used to assemble the obtained clean reads. After assembly, a total of 55,894 unigenes were obtained, with a length range of 201–22,885 bp (mean, 1202 bp; N50, 1964 bp), and a mean GC content of 41.09%. The assembled transcript length distribution is illustrated in Fig. 1.
Unigene function annotation.
Using BLAST X, the assembled unigenes were compared to the protein databases, Nr, SWISS-PROT, KEGG, and KOG (Evalue < 0.00001), to obtain protein function annotation. A total of 33987 unigenes were annotated, of which 33921 were annotated in the Nr database, accounting for 99.8% of all annotated genes; 20750 were annotated in the SWISS-PROT database, 61.1% of all annotated genes; 13165 were annotated in the KEGG database, 38.7% of all annotated genes; and 17462 were annotated in the KOG database, accounting for 51.4% of all the annotated genes.
Screening for DEGs.
DEGs in resistant and susceptible varieties at different time points were screened using P < 0.01 and FC ≥ 2 as the threshold (Fig. 2A). DEGs in the two cultivars at the three time points were compared using Venn diagrams, demonstrating that 2455, 2267, and 1885 DEGs were produced by the susceptible cultivar at 8, 17 and 23 dpi, respectively, among which 1969, 1657, and 1375 DEGs were unique to each time point (Fig. 2B). In the resistant variety, 2324, 4820, and 3584 DEGs were identified at 8, 17, and 23 dpi, respectively, among which 1612, 3621, and 2450 were unique to each time point (Fig. 2C).
Next, DEGs in resistant and susceptible varieties at 17 dpi, relative to before inoculation, were compared using a Venn diagram of DEGs for gb0d-vs-gbjz17d and kb0d-vs-kbjz17d (representing DEGs at GBjz17d and KBjz17d, respectively, relative to those at GB0d and KB0d). The results showed that 8,678 DEGs were unique to the resistant variety at 17 dpi (Fig. 2D), while 2,340 were unique to the susceptible variety, and 2,330 common to both varieties.
GO enrichment analysis of DEGs.
GO enrichment analysis of DEGs in susceptible mulberry at 17 dpi.
Comparing the 2267 DEGs in the susceptible mulberry cultivar at 17 dpi (GBjz17d) with the GO database (Fig. 3A) indicated that they were enriched for 1129 GO entries. There were 747, 280, and 102 GO items in the three major ontologies of biological process, molecular function, and cell component, respectively. The number of down-regulated genes in each item was significantly higher than that of up-regulated genes. In the biological process ontology, the proportions of genes were higher in the classifications, metabolic process, cellular process, and single-organism process. In the ontology of molecular function, the proportion of genes were higher in the classifications, catalytic activity and protein binding. In the cell component ontology, the proportion of genes was higher in the cell, cell part, membrane, organelle, and membrane part classifications (Fig. 3A).
GO enrichment analysis of DEGs in resistant mulberry at 23 dpi.
Comparing the 3584 DEGs in the resistant mulberry cultivar at 23 dpi (KBjz23d) with the GO database (Fig. 3B) demonstrated that they were enriched in 1309 GO entries, including 898, 292, and 119 GO entries in the three major ontologies of biological process, molecular function, and cell composition, respectively. The number of up-regulated genes in each entry was significantly higher than that of down-regulated genes. In the three ontologies, the GO classifications of enriched DEGs for KBjz23d were the same as those for GBjz17d; however, the enriched GO entries differed from those for GBjz17d; in the biological process ontology, KBjz23d DEGs were significantly enriched for entries related to the regulation of macromolecular metabolism, while GBjz17d DEGs were significantly enriched for entries related to plant cell wall metabolism.
GO enrichment analysis of DEGs between resistant and susceptible varieties.
Comparing samples from resistant and susceptible varieties at 17 dpi, there were 2340 unique DEGs (referred to as KBvsGB-GBjz17d DEGs) in susceptible varieties. Comparison of these DEGs with the GO database showed they were mainly enriched for 1146 GO items (Fig. 4A), including 776, 258, and 112 in the ontologies biological process, molecular function, and cell composition, respectively. The number of down-regulated genes was significantly higher than that of up-regulated genes. In the biological process ontology, the proportion of genes was higher in the classifications metabolic process, cellular process, and single organization process, while for the molecular function ontology, the proportion of genes was higher in the classifications catalytic activity and protein binding, and in the cell component ontology, the percentage of genes in the classification cells and cell parts was higher (Fig. 4A).
Next, the 8678 unique DEGs in resistant varieties (KBvsGB-KBjz17d DEGs) were compared with the GO database. The results suggested (Fig. 4B) that KBvsGB-KBjz17d DEGs were mainly enriched in 1764 GO entries, among which there were 1179, 422, and 163 GO entries in the three major ontologies, biological process, molecular function, and cell composition, respectively. The number of down-regulated genes in each entry was significantly higher than that of up-regulated genes. In the three ontologies, the GO classifications enriched for KBvsGB-KBjz17d DEGs was the same as that of KBvsGB-GBjz17d (Fig. 4), while the enriched GO entries differed from those of KBvsGB-GBjz17d. In the cell component ontology, KBvsGB-KBjz17d DEGs were significantly enriched for items including ribosomes and cell-matrix (Fig. 4C), while those of KBvsGB-GBjz17d were significantly enriched for items related to photosynthesis, such as thylakoid membranes (Fig. 4D).
Pathway enrichment analysis of DEGs.
Pathway analysis is helpful for understanding gene biological functions, and for identifying key signal transduction and biochemical metabolism pathways involving DEGs.
Pathway enrichment analysis of DEGs in susceptible mulberry at 17 dpi.
Pathway annotation of DEGs in susceptible mulberry at 17 dpi (GBjz17d) using the KEGG database showed that they were enriched in 113 pathways, with five significantly enriched pathways (Q ≤ 0.05) (Fig. 5A), as follows: biosynthesis of secondary metabolites, ubiquinone and other terpenoid-quinone biosynthesis, flavonoid biosynthesis, phenylpropanoid biosynthesis, and metabolic pathways. Among them, the pathway including the highest number of enriched genes, was the metabolic pathway, which included 145 DEGs. The enrichment factor (rich factor) value of the flavonoid synthesis pathway was the highest, indicating that the pathway was the most enriched.
Pathway enrichment analysis of DEGs in resistant mulberry at 23 dpi.
DEGs from resistant mulberry at 23 dpi (KBjz23d) were annotated with pathways using the KEGG database. The results showed that DEGs were enriched in 109 pathways, six of which were significantly enriched (Q ≤ 0.05) (Fig. 5B): biosynthesis of secondary metabolites; glucosinolate biosynthesis; cutin, suberine, and wax biosynthesis; alpha-linolenic acid metabolism; phenylpropanoid biosynthesis; and flavonoid biosynthesis. Among these, the pathway enriched for the most genes was secondary metabolite synthesis, including 107 DEGs. The glucosinolate synthesis and cutin, suberine, and wax biosynthesis pathways showed the highest level of enrichment.
Pathway enrichment analysis of DEGs between resistant and susceptible varieties.
DEGs with expression differences between resistant mulberry before inoculation and susceptible mulberry at 17 dpi (KBvsGB-GBjz17d) were annotated for pathway enrichment using the KEGG database, demonstrating that DEGs were enriched for 111 pathways, of which two were significantly enriched (Q ≤ 0.05) (Fig. 6A): ubiquinone and other terpenoid-quinone biosynthesis and cutin, suberine and wax biosynthesis. Cutin, suberine and wax biosynthesis was the most strongly enriched.
Annotation of DEGs with expression differences between resistant mulberry before inoculation and at 17 dpi (KBvsGB-KBjz17d) by KEGG database analysis showed that they were enriched for 126 pathways, with seven significantly enriched pathways (Q ≤ 0.05) (Fig. 6B), as follows: starch and sucrose metabolism; pentose and glucuronate conversion interconversions; ABC transporters; phenylpropanoid biosynthesis; plant-pathogen interaction; arginine and proline metabolism; and cyanoamino acid metabolism. Among these pathways’ ABC transporters was the most enriched.
DEG expression trends.
The expression trends of DEGs between susceptible mulberry at 17 dpi (GBjz17d), resistant mulberry at 23 dpi (KBjz23d), resistant mulberry before inoculation and susceptible mulberry at 17 dpi (KBvsGB-GBjz17d), and resistant mulberry before inoculation and susceptible mulberry at 17 dpi (KBvsGB-KBjz17d) were analyzed, resulting in identification of 20 trends. In the trend graphs, colors indicate significantly enriched trends, while those without color are non-significantly enriched trend blocks (Fig. 8–12).
GBjz17d DEG expression trends.
Trend analysis of GBjz17d DEGs identified profile2, profile1, profile18, profile16, profile17, and profile3 (Fig. 7A) as significantly enriched trends, including 452, 312, 210, 219, 196, and 169 genes, respectively (Fig. 7B).
KBjz23d DEG expression trends.
Trend analysis of KBjz23d DEGs identified profile17, profile2, profile19, profile0, profile12, and profile15 (Fig. 8A) as significantly enriched trends, including 627, 260, 349, 333, 260, and 342 genes, respectively (Fig. 8B).
Expression trends of DEGs between resistant and susceptible varieties.
Trend analysis KBvsGB-GBjz17d DEGs (Fig. 9A) revealed profile2, profile17, profile1, profile16, profile18, profile7, and profile0 as significantly enriched and including 461, 347, 244, 307, 227, 136, and 106 genes, respectively (Fig. 9B).
Trend analysis of KBvsGB-KBjz17d DEGs (Fig. 10A) showed that profile2, profile1, profile17, profile18, profile12, profile3, profile7, and profile0 were enriched, including 2003, 953, 1373, 492, 346, 423, 294, and 306 genes, respectively (Fig. 10B).
Co-expression network analysis of key trends.
The core gene, with the highest connectivity in GBjz17d profile1, was unigene 0015188; the co-expression network diagram for this gene is presented in Fig. 11A. Unigene 0015188 encodes β-galactosidase, which is involved in the degradation of pectin in plant primary cell walls. The core gene with the highest connectivity in GBjz17d profile3 was unigene 0022838, and the related co-expression network is presented in Fig. 11B. Unigene 0022838 encodes CASP like protein 1C1, which is related to plant stress resistance. Core genes with high connectivity in the GBjz17dprofile16 were unigenes 0022247 and 0010321; the related co-expression network is shown in Fig. 11C. Unigene 0022247 encodes 3-hydroxy-3-methylglutaryl coenzyme A reductase, which is a key enzyme for antitoxin and steroid based synthesis of sesquiterpenoids, while unigene 0010321 may encode a leucine rich repeat receptor protein kinase (LRR-RLK), which is involved in plant hormone signal transduction, recognition, and transmission between plants and pathogens. The core gene with high connectivity in GBjz17d profile18 was unigene 0021656, and the related co-expression network is detailed in Fig. 11D. Unigene 0021656 encodes calmodulin dependent protein kinase type I, and KEGG annotation showed that this gene is part of the plant-pathogen interaction pathway.
The core gene with highest connectivity in kbjz23d profile0 was unigene 0014848; the related co-expression network diagram is shown in Fig. 11E. Unigene 0014848 encodes ubiquitin carboxy terminal hydrolase 12, which is an effector protein that can inhibit the resistance of plants to root-knot nematodes and promote the formation of nematode feeding points. In kbjz23d profile12, the core gene with highest connectivity was unigene 0015083; the related co-expression network is shown in Fig. 11F. Unigene 0015083 encodes a 16-o-methyltransferase ROMT protein involved in plant alkaloid synthesis. The core genes with higher connectivity in the kbjz23d profile15 were unigenes 0073558 and 0073272, and the related co-expression network is shown in Fig. 11G. Unigene 0073558 encodes ubiquitin carboxyl terminal hydrolase 13, while unigene 0073272 encodes the srg1 protein, which belongs to the zinc finger transcription factor family. The expression of srg1 is related to the accumulation of nitric oxide during plant immune responses. Further, GO biological process annotation showed that this gene is related to flavonoid synthesis. The core genes with higher connectivity in kbjz23d profile19 were unigenes 0002889 and 0068389, the related co-expression network diagram is shown in Fig. 11H. Unigene 0002889 encodes protein phosphatase 2C, which can regulate abscisic acid (ABA) signal transduction and unigene 00683889 encodes a protein kinase. Cipk14, which plays a role in stress responses; GO biological annotation indicated that this gene is related to cell process regulation.
The core gene with higher connectivity in kbvsgb-kbjz17 dprofile3 was unigene 0004006, and the related co-expression network is shown in Fig. 11I. Unigene 0004006 encodes an ERF transcription factor, which is involved in signal transmission using salicylic acid (SA), ethylene, jasmonic acid (JA), and other factors, and plays an important role in plant stress response. In kbvsgb-kbjz17dp profile12, core genes with higher connectivity were unigenes 0000628 and 0015083, and the related co-expression network diagram is presented in Fig. 11J. Unigene 0000628 encodes the MYB transcription factor family protein, AIM1, which contributes to stress-related gene expression during ABA signal transmission and positively regulates plant defense responses to pathogens. Unigene 0015083 encodes a rumt glycyrrhizin-16-o-methyltransferase gene, associated with plant alkaloid synthesis. The core genes with higher connectivity in kbvsgb-gbjz17d profile16 were unigenes 0052595 and 0004809, and the related co-expression network is shown in Fig. 11K. Unigene 0052595 encodes the autophagy-related factor, Atg8, which is under threat in plants. Unigene 0004809 encodes cell wall-related receptor kinase 22, which belongs to the RLK receptor kinase family involved in plant defense.
RT-qPCR validation of DEGs.
To verify the accuracy of the transcriptome sequencing results, 21 DEGs were analyzed by RT-qPCR using cDNA obtained by reverse transcription of the RNA used for sequencing as the template. Relative gene expression levels were calculated using the 2−ΔΔCt method and GraphPad Prism 8.0 used to visualize the results (Fig. 12). Although the magnitude of differences in expression detected by RT-qPCR were not identical to those of DEGs detected by transcriptome sequencing, the direction of the change of DEG expression was consistent between the two approaches, indicating that the results of transcriptome sequencing were highly reliable.