Transcriptomic analysis of gene expression of Verticillium dahliae upon treatment of the cotton root exudates

Background Cotton Verticillium wilt is one of the most devastating diseases for cotton production in the world. Although this diseases have been widely studied at the molecular level from pathogens, the molecular basis of V. dahliae interacted with cotton has not been well examined. Results In this study, RNA-seq analysis was carried out on V. dahliae samples cultured by different root exudates from three cotton cultivars (a susceptible upland cotton cultivar, a tolerant upland cotton cultivar and a resistant island cotton cultivar) and water for 0 h, 6 h, 12 h, 24 h and 48 h. Statistical analysis of differentially expressed genes revealed that V. dahliae responded to all kinds of root exudates but more strongly to susceptible cultivar than to tolerant and resistant cultivars. Go analysis indicated that ‘hydrolase activity, hydrolyzing O-glycosyl compounds’ related genes were highly enriched in V. dahliae cultured by root exudates from susceptible cotton at early stage of interaction, suggesting genes related to this term were closely related to the pathogenicity of V. dahliae. Additionally, ‘transmembrane transport’, ‘coenzyme binding’, ‘NADP binding’, ‘cofactor binding’, ‘oxidoreductase activity’, ‘flavin adenine dinucleotide binding’, ‘extracellular region’ were commonly enriched in V. dahliae cultured by all kinds of root exudates at early stage of interaction (6 h and 12 h), suggesting that genes related to these terms were required for the initial steps of the roots infections. Conclusions Based on the GO analysis results, the early stage of interaction (6 h and 12 h) were considered as the critical stage of V. dahliae-cotton interaction. Comparative transcriptomic analysis detected that 31 candidate genes response to root exudates from cotton cultivars with different level of V. dahliae resistance, 68 response to only susceptible cotton cultivar, and 26 genes required for development of V. dahliae. Collectively, these expression data have advanced our understanding of key molecular events in the V. dahliae interacted with cotton, and provided a framework for further functional studies of candidate genes to develop better control strategies for the cotton wilt disease.


Background
Verticillium dahliae (V. dahliae), a fungal pathogen causing Verticillium wilt, is extremely persistent in the soil and has a broad host range [1,2]. Microsclerotia of V. dahliae overcome the mycostatic activity of the soil and germinate towards roots in the presence of root exudates [3]. The hyphae enter host plants by formation of an infection structure, known as hyphopodium, to develop a penetration peg to pierce root epidermal cells [4]. They enter and clog the xylem vessels, resulting in leaf curl, necrosis, defoliation, vascular tissue wilt, and discoloration [5]. During its life cycle, cotton is continuously threatened by V. dahliae. More than half of the cotton fields in China are affected by V. dahliae and can lead to 30-50% reduction in yield, and even totally wipe out the crop. Verticillium wilt is one of the most severe cotton diseases not only in China but also in other countries. Outbreak of the disease causes substantial economic loss due to significant reduction in fiber yield and quality.
To combat the challenge of V. dahliae, resistance cotton has evolved multiple layers of defense mechanisms, including tissue composition, physiological and biochemical resistance, during the long time period of coexistence and arm race [6][7][8][9][10]. In recent years, with the application of genomics, transcriptomics and proteomics, great progress has been made in understanding the molecular mechanism underlying cotton's resistance against V. dahliae, and a number of genes related to V. dahliae resistance have been identified [11][12][13][14][15][16]. On the other hand, in view of the co-evolving relationship between cotton and V. dahliae, it is also of vital importance to study the molecular mechanisms determining the pathogenicity of V. dahliae. With the completion of genome sequencing of V. dahliae and the development of bioinformatics tools, genomic and transcriptomic sequence information of V. dahliae provide us opportunity for better understanding the pathogenicity of V. dahliae. Analyses of V. dahliae transcriptomes during microsclerotia formation and early infection stage have given us a snapshot of the genes important for development, microsclerotia formation and infection of V. dahliae [17][18][19][20]. For instance, VdPKAC1, VMK1, VdMsb, VdGARP1, VDH1, Vayg1 and VGB were found to be involved in the microsclerotia formation and pathogenic process of V. dahliae [3,[21][22][23][24][25][26]; VdSNF1 and VdSSP1 are related to cell wall degradation [27,28]; VdNEP, VdpevD1, VdNLP1 and VdNLP2 encode effector proteins are involved in the pathogenic reaction [29][30][31][32]; VdFTF1, Vta2 and VdSge1 encode transcriptional factors regulating pathogenic genes [33][34][35]. However, due to the complexity of the pathogenic molecular mechanism of V. dahliae, we still know little about the role of these genes in the interaction between V. dahliae and cotton.
Successful pathogens must be able to recognize and overcome host-plant defense responses [36]. V. dahliae invades cotton through the root system [4,37], therefore, the biological effect of the root exudates is expected to be crucial for successful infection of V. dahliae. Not surprisingly, root exudates have been found to be closely related to plant resistance [38,39]. The root exudates of cotton are rich in amino acids and sugars. Compared with the root exudates from the susceptible cotton cultivars, the root exudates from resistant cotton cultivars lacked aspartic acid, threonine, glutamic acid, alanine, isoleucine, leucine, phenylalanine, lysine and proline, but contained arginine that was absent in the susceptible cottons. No significant difference of saccharide was found in the root exudates between the susceptible and resistant cultivars, but the root exudates of the susceptible cultivars had a much higher concentrations of glucose, fructose and sucrose than that of the resistant ones [40]. Root exudates from the resistant and susceptible cottons inhibited and promoted the growth of V. dahliae, respectively [40][41][42]. However, we know nothing about the molecular basis behind this observation.
In this study, we investigated the effects of root exudates from cotton cultivars susceptible, tolerant or resistant to V. dahliae on the development of the pathogen and performed a time course expression analysis of V. dahliae genes using RNA-seq to (1) compare transcriptomic profiles of V. dahliae in response to root exudates from cottons with different level of V. dahliae resistance, (2) identify biological processes in V. dahliae affected by different root exudates based on analysis of Gene Ontology (GO) terms of the differentially expressed genes, and (3) identify genes involved in the initial steps of roots infection and likely in pathogenesis of V. dahliae. We expect that identification of pathogenic genes in V. dahliae would provide us clues to develop novel strategies for breeding novel cotton germplasm resistant to V. dahliae and/or effective crop management schemes to minimize the infection of V. dahliae.

Methods
Cotton cultivars and V. dahliae strain Two Upland cotton (G. hirsutum L.) cultivars Xinluzao 8 (X) and Zhongzhimian 2 (Z), and one Sea island (G. barbadense L.) cultivar Hai7124 (H) used in this study were collected from the Institute of Cotton Research of Chinese Academy of Agricultural Sciences (Anyang, China) and Shihezi Academy of Agricultural Sciences (Shihezi, China). The 3 cotton cultivars were authorized for only scientific research purpose, and were deposited in the original institutes and College of Agriculture in Shihezi University. The highly virulent V. dahliae strain, V991, was provided and confirmed by the Institute of Cotton Research of Chinese Academy of Agricultural Sciences (Anyang, China). The growth conditions of the cotton cultivars, the preparation of V.dahliae spore suspensions for infection assays and determination of Disease Index after inoculation were described previously [43,44].

Collection of root exudates
Xinluzao 8, Zhongzhimian 2 and Hai7124 are susceptible, tolerant and resistant to V. dahliae, respectively. Cotton seeds were surface sterilized by immersion in 1% (w/v) NaClO and rinsed three times with sterile distilled water. After germination in petri dish, the seeds were sown in sand that were treated by soaking in dilute suphuric acid and sterilized by high temperature. For each cultivar 18 germinated seeds were evenly planted in 2 pots and were grown in a greenhouse with a photoperiod of 16 h light/8 h darkness at 28°C. The cotton seedlings were fed with Hoagland nutrient solution every 3 days (3d). After 45d, the plants were removed from sand, and the sand was immersed with 2 L distilled water to sufficiently dissolve root exudates. The water solution was then filtered with a bacterial filter (0.22 μm in diameter) and concentrated to 0.5 L in a freeze dryer.
V. dahliae strain culture V. dahliae strain, V991, was maintained in 20% glycerol at − 80°C at the Key Laboratory of Oasis Eco-agriculture in Shihezi University. The stored conidia of V991 were incubated on a potato-dextrose agar plate for 1 week and then inoculated into Czapek broth for 5d at 25°C 180 rpm under dark donditions. The fresh conidia and spores were then collected to be used in the root exudate treatment experiments. For each cultivar, 0.5 g of V991 conidia and spores were suspended in 5 mL of root exudates. After cultured for 6, 12, 24 or 48 h at 25°C 220 rpm in 10 mL centrifugal tubes, V991 conidia and spores (Vd-X-6, Vd-X-12, Vd-X-24, Vd-X-48, Vd-Z-6, Vd-Z-12, Vd-Z-24, Vd-Z-48, Vd-H-6, Vd-H-12, Vd-H-24 and Vd-H-48) were collected for RNA extraction. The same amount of V991 conidia and spores suspended in water and cultured for 0, 6, 12, 24 or 48 h were done in parallel (Vd-0, Vd-W-6, Vd-W-12, Vd-W-24 and Vd-W-48). Each time point had two biological replicates. In total, 34 samples were collected and used in RNA-seq.

RNA extraction
Total RNA of V. dahliae was isolated using the RNA simple total RNA kit (Tiagen, Beijing, China) according to the manufacturer's protocol. All RNA samples were treated with RNase-free DNase I. Degradation and contamination of RNA were assessed by using agarose gel electrophoresis. The RNA purity and integrity were determined by a NanoDrop® 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The RNA concentration was measured by a Qubit® 2.0 Fluorometer (Thermo Scientific, Wilmington, DE, USA). High quality RNA samples were chosen for RNA-Seq analyses.

RNA-Seq library construction and sequencing
RNA-Seq library preparation and sequencing were performed at Novogene Bioinformatics Technology Co., Ltd. (Beijing, China) using the standard Illumina protocols. Briefly, mRNAs were enriched from 1.5 μg total RNA by using magnetic beads with Oligo (dT), and then fragmented by adding fragmentation buffer. The short fragments were used as templates to synthesize the first stranded cDNAs with random hexamers. Doublestranded cDNAs were then synthesized by using DNA Polymerase I and RNase H and purified with AMPure XP beads. The purified double-stranded cDNAs was then end repaired, added A tail and ligated with sequencing adapters. The products were enriched with PCR to create the final cDNA libraries. Finally, the library was sequenced on the Illumina Hiseq™ 4000 platform (Illumina, San Diego, CA, USA, 2010).

RNA-Seq data analysis and identification of differentially expressed genes
Raw reads were pre-processed by removing low quality sequences and adaptor using Trimmomatic [45]. The Q30 values, GC content, and sequence duplication levels were calculated for the clean data. All downstream analysis used the clean data with high quality. The resulting high-quality clean reads were then aligned to V. dahliae sequence from genome database (http//www.broadinstitute.org/annotation/genome/Verticillium dahliae/ Blast.html) using the HISAT software [46]. Following alignment, raw read counts for each V. dahliae gene were generated and normalized to FPKM (fragments per kilobase of exon model per million mapped fragments) [47]. The expression level of each gene was analyzed using the union model implemented in the HTSeq software [48]. Differentially expressed genes (DEGs) were identified by using the DEGseq software with the following criteria: a fold change> 2.0 and an adjusted p value< 0.05 [49]. Gene ontology (GO) term enrichment analysis of DEGs was performed based on the Wallenius noncentral hyper-geometric distribution using the GOseq software [50].

qRT-PCR confirmation of differentially expressed genes
Total RNA from V. dahliae was isolated as mentioned above. One microgram of total RNA was used for firststrand cDNA synthesis with the M-MLV reverse transcriptase (TaKaRa, Dalian) according to the manufacturer's instructions. The cDNAs were then used as templates for quantitative real-time PCR (qRT-PCR) experiments. The gene specific primers used in qRT-PCR are listed in Table 1, and the V. dahliae tubulin gene was used as an internal control. The qRT-PCR assays were performed with SYBR Premix Ex Taq (TaKaRa) on a LightCycler 480 system (Roche, USA). All reactions were measured in triplicate. The relative expression ratio of each gene was calculated from the cycle threshold (CT) values using the 2 -ΔΔCT method.

Identification of cotton resistance to V. dahliae infection
In this study, three cotton cultivars with different level of V. dahliae resistance were selected for collection of root exudates. As can be seen from the Fig. 1, severe leaf wilt disease symptoms and premature defoliation were  (Fig. 1d). According to our results about identification of cotton resistance to V. dahliae and previous reports [43,51], Xinluzao 8, Zhongzhimian 2 and Hai7124 were used as cultivars of susceptible, tolerant and resistant to V. dahliae, respectively.

RNA-seq and transcriptome profiles of V. dahliae
To explore the transcriptomic profiling of V991 interacting with root exudates from cotton cultivars with different level of V. dahliae resistance, we generated a total of 34 RNA-seq datasets, 24 from V. dahliae treated by cotton root exudates (Vd-X-6, Vd-X-12, Vd-X-24, Vd-X-48, Vd-Z-6, Vd-Z-12, Vd-Z-24, Vd-Z-48, Vd-H-6, Vd-H-12, Vd-H-24 and Vd-H-48, each with two replicates), 8 from V. dahliae treated by water (Vd-W-6, Vd-W-12, Vd-W-24 and Vd-W-48, each with two replicates) and 2 from untreated V. dahliae, i.e. Vd-0. An overview of the sequencing results is outlined in Table 2. After discarding the low-quality reads, the total number of clean reads per library ranged from 13 to 22 million, and clean bases ranged from 1.97 to 3.22 Gb. Between 11,657,068 and 19,529,825 of these reads were uniquely mapped to the V. dahliae reference genome. The genic distribution of the uniquely mapped reads indicated that most reads (>88.2%) were mapped to exons, and the others were distributed between introns (0.2-0.3%) and intergenic regions (6.7-11.6%) (Additional file 3: Table S1). The Pearson's correlation coefficients (R 2 ) of FPKM distribution between the two biological replicates for each sample were high in each treatment (R 2 = 0.945-0.987, p<0.001), indicating a good level of reproducibility of the RNA-seq data (Additional file 1: Figure S1). The RNA-seq results were also confirmed to be reliable by qRT-PCR using 8 randomly selected genes ( Table 1, Fig. 2) (Additional file 2: Figure S2). For example, the expression levels of these genes peaked at 6 h in Vd-X, but showed no obvious change in Vd-H and Vd-W.
Based on hierarchical clustering using the FPKM values of all genes, it was found that the 17 samples were classified into two groups (Fig. 3). Group I contained all the Vd-6 (Vd-X-6, Vd-Z-6, Vd-H-6 and Vd-W-6) and Vd-12 (Vd-X-12, Vd-Z-12, Vd-H-12 and Vd-W-12) samples as well as Vd-H-24 and Vd-W-24. The expression profiles of these 10 samples were close to that of Vd-0 (CK), which was also clustered in group I. Group II contained all the four Vd-48 (Vd-X-48, Vd-Z-48, Vd-H-48 and Vd-W-48) samples and two Vd-24 (Vd-X-24 and Vd-Z-24) samples. The clustering tree indicated that the gene expression patterns of the two early time points (Vd-6 and Vd-12) were very similar but clearly different from that of the latest time point (Vd-48). The four Vd-24 samples were clustered into the two groups, but were distinct from other samples in the same group by forming a sub-group, suggesting that 24 h could be a transition point regarding the effect of root exudates on the growth of V. dahliae.

Identification of differentially expressed genes (DEGs)
DEGs would offer insights into the metabolic and regulatory changes in V. dahliae when interacting with root exudates from cottons with different V. dahliae resistance, we thus identified DEGs (p<0.05, fold change >2.0) in each interaction using Vd-0 (CK) as a control. Regarding the treatments (root exudates or water), the largest number of DEGs was found in Vd-X vs CK (4602), followed by Vd-Z vs CK (3896), Vd-H vs CK (3227), and Vd-W vs CK (2392) ( Table 3), suggesting that V. dahliae responded to all kind treatments, but responded more strongly to root exudates from the susceptible cultivar (X) than to those from the tolerant (Z) and resistant cultivars (H). Regarding the effect of treated time, the general trend for Vd-X vs CK, Vd-H vs CK and Vd-W vs CK was that the number of DEGs increased with the increased time of treatment, but for Vd-Z vs CK, there were more DEGs at 24 h than other time points. In all three treatments with root exudates, it seemed there were more up-regulated DEGs than down-regulated DEGs at 6 h, but more down-regulated DEGs than upregulated ones at other time points (12 h, 24 h and 48 h) ( Table 3).
To determine the genes of V. dahliae interacted with root exudates, the up-regulated genes in Vd-X, Vd-Z, Vd-H and Vd-W samples in the group I were examined, respectively. By combining up-regulated DEGs in Vd-6 vs CK and Vd-12 vs CK, a total of 339, 302, 327 and 168 DEGs were acquired in Vd-X vs CK, Vd-Z vs CK, Vd-H vs CK and Vd-W vs CK, respectively ( Fig. 4a, b, c, d). These DEGs (339, 302, 327, 168) were combined together to get 631 DEGs (Fig. 4e). Although Vd-H-24 h and Vd-W-24 h were clustered in the groupI, they were analyzed separately because the number of up-regulated genes in Vd-H-24 h (422) and Vd-W-24 h (301) were obviously greater than other samples in group I (Table 3). By combining up-regulated DEGs in Vd-H-24 h vs CK (422) and Vd-W-24 h vs CK (301), a total of 580 DEGs were obtained (Fig. 4f).
The up-regulated genes in Vd-X, Vd-Z, Vd-H and Vd-W samples in the group II were also examined, respectively. By combining up-regulated DEGs in Vd-24 vs CK and Vd-48 vs CK, a total of 1301 and 1283 DEGs were acquired in Vd-X vs CK and Vd-Z vs CK (Fig. 4g, h), respectively. When these DEGs (1301, 1283) were combined together with DEGs in Vd-H-48 vs CK (716) and Vd-W-48 vs CK (367) comparisons, a total of 1652 DEGs were obtained (Fig. 4i).

Gene ontology analyses of DEGs
To further understand the function of these DEGs, we performed gene ontology (GO) analyses to classify the up-regulated genes in group I and group II samples, respectively. For the group I, up-regulated DEGs (631) were mainly enriched in molecular function category ( Fig. 5a; Additional file 4: Table S2). 'hydrolase activity, hydrolyzing O-glycosyl compounds' (p = 1.22E-05), 'hydrolase activity, acting on glycosyl bonds' (p = 2.15E-05) and 'oxidoreductase activity' (p = 0.000309) were the top three significantly enriched terms in the molecular function category. 'transmembrane transport' (p = 3.77E-05), 'carbohydrate metabolic process' (p = 0.001034), 'oxidation-reduction process' (p = 0.001933) were the top three significantly enriched terms in the biological process category. 'extracellular region' (p = 0.000219) is  Table S2) were similar to that of 631 DEGs (Fig. 5a), suggesting that Vd-H-24 and Vd-W-24 were at the same stage of V. dahliae development as the other samples in group I. Therefore, it can be inferred that the response of V. dahliae to island cotton was more prolonged compared with upland cotton. For DEGs (1652) that were up-regulated in the group II, the GO terms changed greatly compared with the group I ( Fig. 5c; Additional file 4: Table S2). These DEGs were mainly enriched in biological process category. 'translation' (p = 1.67E-10), 'peptide biosynthetic process' (p = 4.18E-10) and 'peptide metabolic process' (p = 6.47E-10) were the top three significantly enriched terms in the biological process category. 'structural constituent of ribosome' (p = 3.70E-12) was the most significantly enriched term in molecular function category. 'ribosome' (p = 6.66E-12) and 'ribonucleoprotein complex' (p = 1.81E-06) were the significantly terms enriched in the component category.
It was notable that some genes were related to hydrolase activity, hydrolyzing O-glycosyl compounds and transmembrane transport which have been reported to Fig. 2 The qRT-PCR analyses of the expression of 8 DEGs selected from all DEGs. The 8 DEGs included VDAG_03038 encoding periplasmic trehalase, VDAG_03526 encoding Alpha-glucuronidase, VDAG_04513 encoding hexose transporter protein, VDAG_05015 encoding betagalactosidase, VDAG_07563 encoding sugar transporter STL1, VDAG_08212 encoding lactose permease, VDAG_08286 encoding alpha-glucosides permease MPH2/3, VDAG_09088 encoding MFS transporter. The V. dahliae tubulin gene (VDAG_10074) was used an internal control. All reactions were measured in triplicate. The expression ratio of the gene was calculated from cycle threshold (CT) values using the 2 -ΔΔCT method be closely related to the pathogenicity of fungi, such as cell wall-degrading enzymes, sugar transporter and MFS transporter [52][53][54][55]. This GO terms were significantly enriched in samples of group I, suggesting that these samples were at the critical stage of V. dahliae-cotton interaction (6 h and 12 h). Therefore, V. dahliae samples at 6 h and 12 h were used for further analysis.
We also performed GO analyses to classify the upregulated genes in Vd-X vs CK (1301), Vd-Z vs CK (1283), Vd-H vs CK (716), Vd-W vs CK (367), respectively ( Fig. 7; Additional file 6: Table S4). As expected, the GO enriched terms of the up-regulated genes in Vd-X vs CK (1301), Vd-Z vs CK (1283), Vd-H-48 vs CK (716), Vd-W-48 vs CK (367) were very similar. It was found that 'translation', 'peptide biosynthetic process' Fig. 3 Hierarchical clustering of samples was performed using FPKM values of all genes identified in V. dahliae. The log10 (FPKM+ 1) values were normalized and clustered. Red and blue bands represent high and low gene expression genes, respectively. The color ranges from red to blue, indicating that log10 (FPKM + 1) is from large to small and 'peptide metabolic process' were the top three significantly enriched terms in the biological process category. 'ribosome', 'ribonucleoprotein complex' and 'intracellular non-membrane-bou' were the top three significantly enriched term in the component category. 'structural constituent of ribosome' and 'structural molecule activity' were the significantly terms enriched in molecular function category. No 'transmembrane transport' and 'hydrolase activity, hydrolyzing O-glycosyl compounds' Go enriched terms were found in these samples of group II, again suggesting that 6 h and 12 h were the critical stage of V. dahliae-cotton interaction, while 24 h and 48 h were not.
Genes response to root exudates from different cotton cultivars in V.dahliae GO analyses for the up-regulated DEGs found that transmembrane transport was the most significantly enriched GO term in Vd-X vs CK (339), Vd-Z vs CK (302), Vd-H vs CK (327) comparisons, but not enriched in Vd-W vs CK (168), suggesting that genes related to this term were closely related to the initial steps of the roots infections. Several other GO enriched terms, 'coenzyme binding', 'NADP binding', 'cofactor binding', 'oxidoreductase activity', 'flavin adenine dinucleotide binding', 'extracellular region' were commonly enriched in Vd-X vs CK (339), Vd-Z vs CK (302), Vd-H vs CK (327), suggesting that genes related to these GO terms were also required for the initial steps of the roots infections. Although the main enriched GO terms were similar, the DEGs were quite different in Vd-X vs CK (339), Vd-Z vs CK (302) and Vd-H vs CK (327). Only 57 genes (Fig. 4e) were found to be commonly up-regulated in Vd-X, Vd-Z and Vd-H at the early stages of interaction. The Heatmap of 57 genes indicated that the expression level of these genes were obviously up-regulated in Vd-X, Vd-Z, and Vd-H at one or two time points of cultured, but not obviously up-regulated in Vd-W (Fig. 8a). These genes were considered as potential candidates for involvement in the initial steps of the roots infections. The 57 genes included 31 genes with known functions (Table 4), and 26 genes with unknown functions. Of 31 genes with known functions, it is notable that 7 genes were related to transmembrane transport ( Fig. 8b; Additional file 7: Table S5), including 4 sugar transporter genes (VDAG_09835, VDAG_02051, VDAG_03649, VDAG_09983), 1 pantothenate transporter liz1 gene (VDAG_02269), 1 DUF895 domain membrane protein gene (VDAG_07864) and 1 Inner membrane transport protein yfaV gene (VDAG_00832) ( Table 4). Few genes have been reported to be related to pathogenicity of V. dahliae, such as a gene encoding cyclopentanone 1,2monooxygenase [18], two genes encoding thiamine transporter protein [56,57]. Functional analysis for these candidate genes may be useful for the study of the molecular basis of V. dahliae interacted with cotton.
Genes response to root exudates from susceptible cotton cultivar in V. dahliae GO analyses for the up-regulated DEGs found that 'hydrolase activity, hydrolyzing O-glycosyl compounds' was the most significantly enriched term in molecular function category in Vd-X (339) (p = 8.78E-05) (Fig. 6a), but not in Vd-Z (302), Vd-H (327), Vd-W (168) (Fig. 6b, c, d) suggesting that genes related to this term would be contribute to the pathogenesis of V. dahliae. A total of 20 genes related to this term were found in Vd-X (339), including 16 genes (1-16) reported to be related to cell wall degradation (Table 5) [58]. A total of 121 DEGs unique to Vd-X (Fig. 4e) whose expression were up-regulated only in root exudates from susceptible cotton cultivar (X) were thought to be the candidate genes related to pathogenesis of V. dahliae. The Heatmap of 121 genes indicated that the expression level of these genes were obviously up-regulated in Vd-X, and only few genes were also up-regulated in Vd-Z, Vd-H, and Vd-W at one or two time points of cultured (Fig. 9a). The 121 DEGs included 68 genes with known functions (Table 6), 57 genes with unknown functions. Of 68 DEGs with known functions, it is notable that 9 genes related to hydrolase activity, hydrolyzing O-glycosyl compounds ( Fig. 9b; Additional file 8: Table S6) encode cell walldegrading proteins, including endo-1,4-beta-xylanase (VDAG_03790, VDAG_06165), xylosidase/arabinosidase (VDAG_01866), mixed-linked glucanase (VDAG_07983), glucanase (VDAG_09516), trehalase (VDAG_03038), Alpha-glucosidase (VDAG_01555), Alpha-glucuronidase (VDAG_03526), Alpha-N-arabinofuranosidase (VDAG_ 03553), 13 genes were related to transmembrane transport, including 6 sugar transporter genes (VDAG_07141, Additionally, GO analysis of 66 DEGs unique to Vd-Z (Fig. 10a) and 109 DEGs unique to Vd-H ( Fig. 10b; Additional file 9: Table S7) did not find hydrolase activity, hydrolyzing O-glycosyl compounds and transmembrane transport enriched GO terms, suggesting that the number of DEGs related to hydrolase activity hydrolyzing O-glycosyl compounds and transmembrane transport in Vd-X vs CK (339) were higher than that in Vd-H vs CK (327) and Vd-Z vs CK (302) and these genes may be related to pathogenesis of V. dahliae.

Genes related to development of V. dahliae
A total of 55 genes (Fig. 4e) whose expression were upregulated in Vd-X, Vd-Z, Vd-H and Vd-W were considered to be required for development of V. dahliae. The Heatmap of 55 genes indicated that the expression level of these genes were obviously up-regulated in Vd-X, Vd-Z, Vd-H and Vd-W at one or two time points of cultured (Fig. 11a), which was consistent with the Veen diagram results. The 55 genes included 26 genes with  (Table 7), it is notable that several genes were associated with FAD binding and RNA processing ( Fig. 11b; Additional file 10: Table S8), such as VDAG_02063, VDAG_05832, VDAG_09806, VDAG_05829, VDAG_02981. Functional analysis for these candidate genes may be useful for the study of the molecular basis of V. dahliae development.

Discussion
V. dahliae can survive for many years in soil and dead plant tissues, making Verticillium wilt difficult to control, which has been likened to a bottleneck in commercial crop productivity [53,56]. Only limited studies have focused on pathogenicity-related molecular mechanisms in the fungus, In this study, RNA-Seq was firstly used to explore and compare the transcriptomic profiles of V. dahliae after cultured with root exudates from different cotton varieties. Statistical analysis of DEGs in V. dahliae samples vs CK (Vd-0) revealed that V. dahliae responded to all kinds of root exudates but was more responsive to susceptible cultivar than to tolerant and resistant cultivars. GO analysis revealed the enriched GO terms of up-regulated genes in Vd-X vs CK (339), Vd-Z vs CK (302), Vd-H vs CK (327) were similar. However, the up-regulated genes were quite different in these samples, and only 57 up-regulated genes were found to be common in Vd-X vs CK (339), Vd-Z vs CK (302) and Vd-H vs CK (327), suggesting that the molecular mechanism of the response of V. dahliae to different root exudates from three cotton cultivars was different. GO analysis also found that enriched GO terms of upregulated genes in Vd-X (339) and Vd-Z (302) at 6 h and 12 h of cultured were obviously different from that of Vd-X (1031) and   Plant pathogenic fungi can produce a range of cell wall-degrading enzymes to facilitate infection and colonization [59,60], including cellulase, hemicellulase, pectinase, etc. Hydrolytic enzymes, particularly cellulases and pectinases, have been considered to be important for the expression of disease symptoms and pathogenesis of V. dahliae [61,62]. The cell wall-degrading enzymes are virulence factors, such as such as xyloglucan-specific endoglucanase [63], fungal endopolygalacturonases [64], and also function as pathogen-associated molecular patterns (PAMPs). Specifically, the cell wall-degrading enzymes contain carbohydrate-binding modules (CBM), non-catalytic protein domains that are generally associated with carbohydrate hydrolases in fungi, which are known to act as elicitors of the PAMP-triggered immunity (PTI) response in oomycetes [65,66]. In V. dahliae, two Glycoside hydrolase 12 (GH12) proteins, VdEG1 and VdEG3 acted as PAMPs to trigger cell death and PTI independent of their enzymatic activity in Nicotiana benthamiana.
Although cell wall-degrading enzymes have been received to be related to pathogenicity of fugus, but the  direct molecular evidence was not sufficient. In this study, GO analyses for the up-regulated DEGs found that genes related to hydrolase activity, hydrolyzing Oglycosyl compounds was the most significantly enriched term in molecular function category for Vd-X (339), but not in Vd-Z (302), Vd-H (327), Vd-W (168), including 16 cell wall-degrading genes, suggesting these genes would be contribute to the pathogenesis of V. dahliae. Additionally, A total of 121 DEGs unique to Vd-X (339) whose expression were obviously up-regulated after cultured with root exudates from susceptible cotton cultivar, including 9 cell wall-degrading genes. These results provided a proof of the involvement of cell walldegrading genes in the initial steps of the roots infections and likely in pathogenesis. Recently, functional studies of cell wall-degrading related genes by targeted gene knockout have been carried out to obtain mutants deficient in one or more these genes [60,67], but were not conclusive due to the multigene families encoding these enzymes [68]. Therefore, it is important to detect which genes were responsible for the pathogenicity of V. dahliae. In this study, 16 cell wall-degrading related genes were significantly up-regulated in Vd-X at early stage of interaction, which can be used as the target genes for studying V. dahliae pathogenicity by gene knockout. Here some genes were up-regulated in V. dahliae cultured by water, maybe resulted from no nutrient in water. Perhaps the starvation of the fungus may induce expression of genes encoding cell walldegrading enzymes [69]. The adaptation of V. dahliae inside the host plants requires a large number of channel proteins to control the absorption of nutrients across the plasma membrane [56]. Transport proteins are integral transmembrane protein that exist permanently within and span the membrane across which they transport substances. GO analyses found that transmembrane transport term was commonly enriched in Vd-X (339), Vd-Z (302), Vd-H (327), but not enriched in Vd-W (168) at 6 h and 12 h of cultured, suggesting that they were required for the initial steps of the roots infections. Seven genes related to transmembrane transport found to be up-regulated in V. dahliae cultured by different root exudates, and 13 genes related to this term were only up-regulated in V. dahliae cultured by root exudates from susceptible cultivar. The results exhibited that genes related to this term can respond quickly to cotton root exudates, especially to the susceptible cotton, suggesting that genes related to transmembrane transport may be associated with the initial steps of the roots infections and likely in pathogenesis. The content of carbohydrate and amount of amino acids in the root exudates of susceptible cultivar was distinctly more than resistant ones [42]. Thus, V. dahliae can obtain more nutrients to provide its growth   in root exudates from susceptible cotton, which may be responsible for the higher expression of transmembrane transport genes at the early stage of interaction in V. dahliae cultured by root exudates from susceptible cotton. However, few transmembrane transport genes for nutrient acquisition have been identified from V. dahliae, and their involvement in the disease process is unknown.
In short, our study firstly revealed the transcriptomes of V. dahliae cultured with root exudates from different cotton cultivars. Our results provided the clear proof at the molecular level for the association of cell walldegrading and transmembrane transport related genes with pathogenesis of V. dahliae. The results enriched the genomic information on V. dahliae in public databases, and laid a foundation for the evaluation and understanding the molecular mechanisms of V. dahliae interacted with cotton and pathogenicity. The paper provided a framework for further functional studies of candidate genes to develop better control strategies for the cotton wilt disease.

Conclusions
In this study, we present the first comparative transcriptomic profiling analysis of V. dahliae responded to root exudates from a susceptible upland cotton cultivar, a tolerant upland cotton cultivar and a resistant island cotton cultivar. Our study provided a comprehensive examination of the biological processes in V. dahliae affected by different root exudates based on analysis of Gene Ontology (GO) terms of the differentially expressed genes, and described genes that were involved in the initial steps of the roots infections and likely in pathogenesis. Genes related to 'hydrolase activity, hydrolyzing Oglycosyl compounds' highly enriched in V. dahliae cultured by root exudates from susceptible cotton at early stage of interaction may be responsible for the pathogenicity of V. dahliae. Genes related to 'transmembrane   transport' enriched in different root exudates, but not in water may be required for the initial steps of the roots infections. These expression data have advanced our understanding of key molecular events in the V. dahliae interacted with cotton, and provided a framework for further functional studies of candidate genes to develop better control strategies for the cotton wilt disease.