Study of the main physiological indexes in sunflower roots post-inoculation
Two sunflower genotypes, previously classified as resistant (S18) and susceptible (P77) to Sunflower Verticillium wilt (SVW) according to their field behavior, were used in this study (Fig. 1a). To investigate the process of SVW colonization, 3 indexes (Soluble protein, Peroxidase (POD), and Malondialdehyde (MDA)) were measured at 6, 12 and 24 h, and 2, 3, 5 and 10 days post inoculation for the two sunflower genotypes. After infection with V. dahliae, the content of soluble protein initially increases, but then decrease thereafter, and the contents peaked from 6 to 12 h when treatment was prolonged. The physiological index in the resistant genotype was significantly higher than that in the susceptible genotype at 6, 12, and 24 h and 2 d. However, in the resistant genotype, the physiological index was lower than that in the susceptible genotype at 3 and 5 d, and this reversed at 10 d (Fig. 1b). The POD activity in the resistant genotype was higher at 6, 12, and 24 h, and 2, 5 and 10 d compared with the susceptible genotype (Fig. 1b). For MDA measurement, the content was slightly higher in the susceptible genotype than in the resistant genotype at the four initial time points (Additional file 1: Figure S1). Based on the above discussion of the physiological indexes, we concluded that the resistant genotype performed better in resistance to V. dahliae infection. To further identify the resistance-related genes, a RNA-Seq experiment for all time points examined above was proposed in the following study.
Illumina sequencing and de novo assembly
In the present study, we performed transcriptome analysis of 16 samples (Fig. 2b and c) to describe the sunflower root response to SVW, and 509,533,702 reads were sequenced for library establishment. The length of most of the reads were distributed between 125 and 175 bp. Thus, we used the short reads comparison software TMAP [20] to compare reads with the reference gene sequences (allowing two-base mismatches), obtaining a total of 492,889,748 mapped reads (96.73% of the total reads) (Additional file 2: Table S1). Moreover, we obtained a total of 76,011 unigenes, with a mean length of 890 bp. The size distribution showed that 33,287 (43.8%) unigenes ranged in size from 200 to 500 bp, 19,344 (25.4%) unigenes ranged in size from 500 to 1000 bp, and 30.8% unigenes showed sizes greater than 1000 bp.
The expression level for each gene was calculated using the RPKM [21] method (Reads Per kb per Million reads). The related information for each gene (coverage, symbol, description, etc.) was also provided. To gain insight into the functions of the genes with responses to SVW, all unigenes were annotated using WEGO [22] and Blast2GO [23] software and were classified into functional categories. KEGG Pathway significant enrichment analyses were also performed for all genes. Specific expression analysis [24] was also used in the present study (Additional file 1: Figure S2). Scatter plot, principal component and cluster analyses were also used to mine the deep relation of different samples (Fig. 2a, b and c). These results revealed a large difference between the resistant and susceptible genotypes.
Inter-genotypes differences in basal gene expression
The uninoculated resistant (S18) and susceptible (P77) genotypes were further compared to analyze the basal gene expression pattern (Fig. 2a). Differential expression analysis revealed 4937 significantly DEGs based on the False Discovery Rate Value (FDR) < = 0.001 and the absolute value of fold-change (FC) > = 2 (Additional file 3: Table S2). Moreover, among all DEGs, 1057 genes had a much higher expression value for the resistant genotype, with RPKMs ranging from 15 to 2685, and a total of 962 genes had a much higher expression value for the susceptible genotype, with RPKMs ranging from 15 to 1194. In particular, 358 genes were expressed at low level in the susceptible genotype (FC <0.001), but they had a higher expression value in the resistant genotype (FC > = 10). As inter-genotype differences may reveal the mechanisms of disease resistance and susceptibility [25], we also classified these DEGs according to functional categories. Following the Nr annotations, the DEGs were mapped into three GO categories, including 58 sub-categories (Additional file 1: Figure S3). Moreover, KEGG analyses were performed to identify the basal level biological pathways in sunflower roots. All DEGs were enriched into 120 KEGG pathways. Moreover, 16 important pathways were obtained with P values < =0.01 (Additional file 4: Table S3, Additional file 1: Figure S4).
Transcriptional changes in response to V. dahliae inoculation
After the multiple comparisons (Additional file 1: Figure S5), two criteria of DEGs were defined: |FC value| > = 2 with P-value < = 0.01, and FDR < = 0.001 [26]. The results of the differential comparison between control and inoculated samples are shown in the RD and SD datasets (Additional file 5: Table S4-6). The 14 datasets represented the DEGs in response to SVW for each genotype. The comparison results between two inoculated samples for resistant and susceptible genotypes at same time point is shown in the D datasets, representing the DEGs between resistant and susceptible genotypes in response to SVW infection (Additional file 5: Table S4–6).
The dynamic changes with time of the DEGs in RD, SD and D datasets were also investigated (Additional file 1: Figure S6). The results showed that the remarkable changes in the transcriptome profile occurred on day 2. The number of DEGs in the D datasets, including both total and up-regulated DEGs, was markedly more than those in other two datasets at the first three time points. For the RD datasets (total, up or down-regulated DEGs), the numbers of DEGs was much more than that in the SD datasets at the last three time points (3, 5, and 10 d). In contrast, the numbers of DEGs in the SD datasets have relatively fewer changes than the other two datasets at the last three time points (3, 5, and 10 d).
Further, time-common DEGs and time-specific DEGs were examined (Fig. 3a, b, c). The common DEGs between neighboring time points were defined as time-common DEGs. The specific DEGs for their own datasets were defined as time-specific DEGs. The time-common DEGs showed more overlap in the datasets of 2–3 d in RD and 24 h-2 d in SD (Fig. 3I, III). For D datasets, the time-common DEGs were more common in neighboring data than in the RD and SD datasets (Fig. 3I, III, V). And, Fig. 3II, IV and VI show the trend of time-specific DEGs.
Genotype-common transcriptional changes in response to V. dahliae inoculation
For the time-common DEGs, the property of the genes between RD and SD datasets was also investigated. A total of 1319 co-expressed DEGs were defined as genotype-common transcriptional DEGs (Additional file 6: Table S7). The functional categories of the largest percentages of the common genes were related to ‘metabolic process’ (17.3%) and ‘cellular process’ (16.5%), followed by ‘disease resistance’ (18.8%) (Fig. 4a). Among these disease resistance categories, the cluster ‘response to stress’ (234, 8.1%) represented the largest group, followed by ‘oxidation reduction’ (115, 3.9%); ‘secondary metabolic process’ (82, 2.8%); and ‘immune response’ (42, 1.4%). Others disease resistance categories, such as ‘hormone metabolic process’ , ‘positive regulation of immune system process’ , ‘regulation of response to stimulus’ and ‘cell wall biogenesis’ , were also studied (Fig. 4c).
The 5956 time-common DEGs in the D datasets were also co-expressed in both genotypes and were also observed as DEGs with genotype-common transcriptional pattern (Additional file 6: Table S7). The functional categories of the largest percentages of the common genes were related to ‘metabolic process’ (17.9%) and ‘cellular process’ (17.9%), followed by ‘disease resistance’ (16.6%), which included the categories: ‘response to stress’ (677, 6.5%); ‘oxidation reduction’ (292, 2.8%); ‘secondary metabolic process’ (178, 1.7%); and ‘immune response’ (142, 1.4%). Other disease resistance categories were also observed following detailed analysis (Fig. 4b, c).
After analyzing the genotype-common transcriptional pattern, 1231 genes related to SVW resistance were observed in the RD-VS-SD and D datasets (Additional file 6: Table S7, Additional file 7: Table S8). The genes encoding to peroxidase (POD), glutathione peroxidase, aquaporin PIP, chitinase, L-ascorbate oxidase, and LRR receptors were identified using GO enrichment (Additional file 7: Table S8). For example, unigene_26949, which encoded peroxidase, showed higher FCs at different time points in the resistant genotype. The inverse result was observed for gene CL7892, which encoded L-ascorbate oxidase (Fig. 4d). When we focused on the average fold-change of each DEG for each functional category in genotype-common DEGs, most of the up-regulated DEGs in the resistant genotype showed a higher average FC than the susceptible genotype, and several down-regulated DEGs in the resistant genotype displayed lower average FCs than the susceptible genotype (Additional file 7: Table S8).
Genotype-specific transcriptional changes in response to V. dahliae inoculation
Certain time-common DEGs were unique to the resistant (S18) or susceptible (P77) genotypes and were defined as genotype-specific transcriptional DEGs. A total of 4112 and 3007 DEGs were observed as genotype-specific transcriptional DEGs in resistant and susceptible genotypes, respectively (Additional file 8: Table S9). The functional enrichment of DEGs showed that the largest percentages of the genes were related to ‘metabolic process’ (17.4 and 17.8%, respectively) and ‘cellular process’ (17.6 and 16.7%, respectively), followed by genes involved in ‘disease resistance’ (17.7 and 18.2%, respectively) (Additional file 1: Figure S7A and B). Among these disease resistance categories in the resistant and susceptible genotypes, most of the DEGs were involved in the cluster ‘response to stress’ (621 and 415, respectively); followed by ‘oxidation reduction’ (274 and 231, respectively); secondary metabolic process’ (177 and 142, respectively); and ‘immune response’ (135 and 74, respectively, Additional file 1: Figure S7C). We concluded that the resistant genotype expressed more disease-resistance genes than the susceptible genotype (Additional file 1: Figure S7).
Analysis of the genotype-specific transcriptional pattern revealed 759 and 511 genes directly related to SVW resistance in the resistant and susceptible genotypes, respectively (Additional file 9: Table S10). Transcripts related to peroxidase (POD), L-ascorbate oxidase, aquaporin PIP, LRR receptor-like serine and the proto-oncogene protein myb were identified using GO enrichment. A total of 22 genotype-specific transcriptional DEGs encoding peroxidase (POD) were identified in the resistant genotype, while only 13 DEGs were identified in the susceptible genotype. Remarkably, for up-regulated genotype-specific DEGs in the resistant genotype, a higher average FC was observed in the resistant genotype than in the susceptible genotype (Fig. 4e, left part of red line). The down-regulated genotype-specific DEGs in the resistant genotype showed a lower average FC (Fig. 4e, left part of red line). Similar results were observed in genotype-specific transcriptional DEGs encoding L-ascorbate oxidase, aquaporin PIP, LRR receptor-like serine kinase and the myb proto-oncogene protein (Additional file 10: Table S11). However, no similar pattern was observed in the susceptible genotype (Fig. 4e, right part of red line). This conclusion was consistent with the results of the analysis of POD physiological indexes in sunflower roots post inoculation (Fig. 1).
Validation of RNA-Seq data by quantitative real-time PCR
To validate the RNA-Seq expression profiles of DEGs, the quantitative real-time PCR (qPCR) of three independent replicatates was performed. A total of 18 DEGs were randomly selected from the comparison between the uninoculated resistant and susceptible genotypes to validate the RNA-Seq expression profiles, and 14 DEGs detected using qPCR were consistent with the RNA-Seq data (Additional file 1: Figure S8). Further, we randomly selected 10 DEGs from RD and SD datasets (Additional file 11: Table S12). For RD datasets, 90.0% (63/70) of the qPCR results were consistent with the RNA-Seq data (Fig. 5a). For SD datasets, 84.3% (59/70) of the qPCR results were consistent with the RNA-Seq data (Additional file 1: Figure S9).
Generally, most of the genes assayed in the inoculated resistant genotype showed higher expression levels, and they were further analyzed by RNA-Seq and confirmed using qPCR, consistent with results of studies of the major physiological indexes in sunflower roots post inoculation (Fig. 5b-i and ii).
Identification of novel V. dahliae inoculation-responsive genes
By analyzing the inter-genotype differences in basal gene expression, genotype-common transcriptional changes and genotype-specific transcriptional changes in response to V. dahlia inoculation, showed 2107 DEGs (Additional file 12: Table S13), which were related to SVW resistance. The results of cluster analysis for the DEGs showed that the expression patterns of the DEGs significantly varied in response to V. dahliae in both genotypes at different time points. Generally, similar expression patterns were observed among adjacent time point samples in the same genotype. In both resistant and susceptible genotypes, the samples possessed more up-regulated genes 2 days after inoculation. Interestingly, Some DEGs were up-regulated in the resistant genotype and down-regulated in susceptible genotype, while other DEGs showed the opposite result (Fig. 2c).
Interestingly, in the interaction mechanism of Verticillium wilt and sunflower, for 2107 DEGs, KEGG enrichment analyses were performed to identify the biological pathways in sunflower root. The results showed that 112 unigenes involved in KEGG pathway related to plant-pathogen interaction, which regulate 35 crucial points (Additional file 1: Figure S10); 97 DEGs were involved in KEGG pathway related to plant hormone signal transduction, which regulate 20 crucial proteins (Additional file 1: Figure S11); 53 DEGs were associated with KEGG pathway related to flavonoid biosynthesis, which regulate 7 critical points (Additional file 1: Figure S12); and 7 DEGs were involved in the KEGG pathway related to benzoxazinoid biosynthesis (Additional file 13: Table S14, Additional file 1: Figure S13).
Based on plant-pathogen interaction pathway, we observed that the hyper-sensitive response (HR) was regulated through the reactive oxygen species (ROS) and oxide (NO) signaling pathways. Respiratory burst oxidase homolog (RBOH) and flagellin-sensitive 2 (FLS2) could activated ROS. Cyclic nucleotide-gated ion channel (CNGC) could mediate the NO signaling pathway (Additional file 1: Figure S14A). In the present study, three genes (Unigene21981_All, Unigene20875_All and Unigene11922_All) encoding RBOH were up-regulated (average FC > 0) in both genotypes, based on the comparison between treated and untreated sunflower genotypes. Notably, the expression levels of the genes in the resistant genotype were highly increased by comparing with those in susceptible genotype after pathogen challenge. However, most of the genes encoding FLS2 were down-regulated in both the resistant or susceptible genotypes after inoculation, and the expression levels of the genes in the resistant genotype were highly decreased by comparing with those in the susceptible genotype. Perhaps to compensate for the loss, the expression levels of three genes in the resistance genotype were more increased. Furthermore, the results showed that the genes encoding CNGC were down-regulated or up-regulated in both genotypes after the infection. However, expression levels of some of genes were more enhanced in the resistant genotype. Remarkably, four genes encoding CNGC were up-regulated in the resistant genotype, but were down-regulated in the susceptible genotype (Additional file 13: Table S14). In plant hormone signal transduction pathway, the dissociation between JAZ (jasmonate ZIM domain-containing protein) and the transcription factor MYC2 could induce defense-related genes, and was initiated by jasmonic acid (JA). This dissociation could be activated by COI1 (Additional file 1: Figure S14B). In this study, the genes (e.g., CL200.Contig5_All) related to COI1 were promoted in the resistant genotype, but were enhanced in the susceptible sunflower genotype when compared treated with untreated sunflowers. Four genes encoding JAZ were up-regulated in the susceptible sunflower genotype, but two genes (Unigene15789_All, CL9889.Contig2_All) were down-regulated in the resistant sunflower genotype. For the salicylic acid (SA)-mediated signal transduction pathway, to induce the expression of defense genes, such as the pathogenesis-related gene PR-1 (pathogenesis-related protein 1), NPR1 must activate the activity of transcription factor TGA. In the present study, Unigene28388_All, which encoded NPR1, was up-regulated in the resistant sunflower genotype, but were down-regulated in the susceptible sunflower genotype. Four genes encoding TGA were up-regulated in the resistance sunflower genotype. Notably, the expression levels of the genes in the resistant genotype were remarkably increased when comparing them with those in the susceptible genotype (Additional file 13: Table S14).
Furthermore, we identified 20 ABA (abscisic acid)-related genes, among which 14 genes were up-regulated in the resistant genotype compared with the susceptible genotype, and 3 DEGs showed distinct down-regulation in the resistant genotype compared with the susceptible genotype. In addition, we identified 8 JA-related genes and 6 ET-related genes (Fig. 6). All these genes were well represented in both genotypes and a higher average FC in the resistant genotype than in the susceptible genotype was observed for many of the up-regulated genes, while a lower average FC in the resistant genotype than in the susceptible genotype was observed for many of the down-regulated genes (e.g. the DEGs encoded membrane proteins, NO and WRKYs), confirming the results of genotype-specific transcriptional pattern analysis (Fig. 6, Additional file 12: Table S13). In particular, these down-regulated genes may regulate important signaling pathway. This conclusion was consistent with the research results of Sun et al., revealing that the down-regulation of GhCYP82D might regulate the JA signaling pathway [9].