Large-scale identification of Gossypium hirsutum genes associated with Verticillium dahliae by comparative transcriptomic and reverse genetics analysis

Verticillium wilt is a devastating disease of cotton, which is caused by the soil-borne fungus Verticillium dahliae (V. dahliae). Although previous studies have identified some genes or biological processes involved in the interaction between cotton and V. dahliae, its underlying molecular mechanism remains unclear, especially in G. hirsutum. In the present study, we obtained an overview of transcriptome characteristics of resistant upland cotton (G. hirsutum) after V. dahliae infection at 24 h post-inoculation (hpi) via a high-throughput RNA-sequencing technique. A total of 4,794 differentially expressed genes (DEGs) were identified, including 820 up-regulated genes and 3,974 down-regulated genes. The enrichment analysis showed that several important processes were induced upon V. dahliae infection, such as plant hormone signal transduction, plant-pathogen interaction, phenylpropanoid-related and ubiquitin-mediated signals. Moreover, we investigated some key regulatory gene families involved in the defense response, such as receptor-like protein kinases (RLKs), WRKY transcription factors and cytochrome P450 (CYPs), via virus-induced gene silencing (VIGS). GhSKIP35, a partner of SKP1 protein, was involved in ubiquitin-mediated signal. Over-expression of GhSKIP35 in Arabidopsis improved its tolerance to Verticillium wilt in transgenic plants. Collectively, global transcriptome analysis and functional gene characterization provided significant insights into the molecular mechanisms of G. hirsutum-V. dahliae interaction and offered a number of candidate genes as potential sources for breeding wilt-tolerance in cotton.


Introduction
Cotton (Gossypium spp.) is a widely cultivated plant and has important economic value because of its fiber and oil. In many cotton diseases, Verticillium wilt is one of the most devastating diseases of cotton worldwide. In China, Verticillium wilt affects more than 70% of cotton planting areas, leading to huge economic loss every year. Verticillium wilt is caused by the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 soil-borne fungus Verticillium dahliae Kleb. It is difficult to control this pathogen, due to its long-term survival as microsclerotia in the soil and broad host range. The primary strategies mainly include breeding and cultivation of resistant varieties. G. hirsutum is known as upland cotton, and it produces more than 95% of the annual cotton crop worldwide. However, the underlying genetic and molecular mechanisms of cotton resistance to Verticillium infection remain poorly explored.
Previous studies have mainly focused on the molecular characterization of the defense responses of G. barbadense upon V. dahliae infection. The phenylpropanoid pathway plays a key role during the plant defense response to V. dahliae [1], which has been also confirmed in tomato [2] and pepper [3]. Xu et al. [1] performed RNA-sequencing in island cotton variety (G. barbadense) and identified the central role of lignin metabolism in cotton resistance to V. dahliae. The lignin level is also detected by histochemical analysis, showing that lignin levels are higher in resistant plants compared with susceptible plants. Ethylene (ET) has been shown to have dual roles in resistance to Verticillium wilt. ET biosynthesis and signaling components are activated in response to V. dahliae infection, which can trigger the systemic resistance in plants [4]. However, former study has shown that released ET makes cotton susceptible to Verticillium wilt, resulting in wilting leaves [5]. It has been identified that several pathogenesisrelated (PR) proteins, such as Bet v1 family and major latex protein (MLP), play important roles in defense reaction against Verticillium wilt [4,6,7,8]. These findings indicated that cotton defense response to V. dahliae infection is sophisticated, which might involve different defense pathways.
Several previous studies have characterized transcriptomic changes in the cotton defense response to V. dahliae infection [1,9,10]. Our group has identified 99 differentially expressed genes (DEGs) in resistant G. hirsutum varieties compared with susceptible varieties by suppression subtractive hybridization (SSH) [7]. But, these data are not enough to elucidate the defense mechanism of cotton against V. dahliae. Moreover, a comprehensive understanding of molecular mechanisms of cotton defense response to V. dahliae has not been established due to complexity of cotton genome. Genome-wide transcriptome analysis is a potentially valuable strategy to elucidate the underlying molecular mechanisms of the physiological processes, and it is substantially efficient to identify genes of interest. Moreover, virus-induced gene silencing (VIGS) offers an efficient approach for large-scale functional analysis of individual genes [11], which will elaborate our understanding of cotton defense response to Verticillium wilt.
In the present study, we aimed to deepen our understanding of upland cotton defense mechanisms against V. dahliae. A comprehensive transcriptome analysis of upland cotton under V. dahliae attack was performed via high-throughput RNA-sequencing technology. Furthermore, in combination with VIGS technique, we further characterized DEGs functions against V. dahliae in upland cotton. Our present data provided valuable information on the molecular basis of cotton defense against V. dahliae and also provided some important candidate genes against V. dahliae.

Plant materials and pathogen inoculation
Two varieties of upland cotton (G. hirsutum), including Zhongzhimian KV3 (Verticillium wilt-resistant upland cotton) and 86-1 (Verticillium wilt-susceptible upland cotton) were used in the present study. Seed germination and seedling growth were carried out at day/night temperatures of 26/20˚C with 1/3 MS medium [12] in sterile culture pots under long-day condition (16-h photoperiod) until the second true leaves appeared. Cotton seedlings were inoculated using root dipping method [7]. Inoculation was performed using the highly toxic and defoliant strain of V. dahliae, V991. Conidial production was prepared based on a previously described method [7], and the concentration of the conidial suspension was adjusted to 1.0×10 7 conidia/mL. The root tissues were harvested at 24 h post-inoculation (hpi). While additional seedlings were mock-inoculated appropriately with sterile distilled water as noninoculated controls. More than 10 plants per treatment were collected and then immediately frozen in liquid nitrogen.

Library construction and Illumina sequencing
Total RNA was extracted using a modified cetyltrimethylammonium bromide (CTAB) method. The integrity of RNA was assessed using denaturing agarose gel electrophoresis and spectrophotometrically examined according to the ratio of A260/A280 (Perkin-Elmer, USA). The RNA samples, including samples of Zhongzhimian KV3 infected with V. dahliae strain V991 at 24 h and mock-inoculated samples as non-inoculated controls, were sent to Beijing Genomic Company (BGI, Shenzhen, China). The cDNA library construction and Illuminabase paired-end deep sequencing were performed following the manufacturer's instructions (Illumina). Two biological replicates were used for all RNA-sequencing experiments. Raw reads produced by Illumina HiSeq 2000 were filtered into clean reads by removing reads containing adapter, unknown bases or low quality reads. After screening, clean reads were aligned to the reference sequences with SOAP aligner/SOAP2 [13]. Distribution of reads was determined according to data alignment, and coverage analysis was then performed.

Gene expression and data analysis
The gene expression level was calculated and normalized by using RPKM (reads per kilobase transcriptome per million mapped reads) method [14]. DEGs of replicates between inoculated and mock-inoculated groups were performed using the DESeq R package (1.10.1) based on the negative binomial distribution. The P-values were adjusted using the Benjamini and Yekutieli's method [15] to control the false discovery rate (FDR). Genes with FDR 0.01 and the absolute value of Log 2 Ratio!2 were used as the threshold to identify the significant DEGs [16] ( Audic et al., 1997). Moreover, functions of DEGs were classified by gene ontology (GO) enrichment analysis and pathway assignments in KEGG (Kyoto Encyclopedia of Genes and Genomes) database [17,18].

qRT-PCR analysis
Total RNA was isolated as previously described. Purified RNA was treated with DNase I (TaKaRa, Japan), then reversely transcribed into cDNA using M-MLV Reverse Transcriptase (Promega, USA) and quantified with spectrophotometry (Perkin-Elmer, USA). Quantitative amplification was performed in triplicate using SuperReal PreMix Plus (Tiangen, China) on an ABI 7500 Real Time PCR System. Briefly, after an initial denaturation step at 95˚C for 15 min, the amplification was carried out with 40 cycles at a denaturating temperature of 95˚C for 10 sec, and an annealing temperature of 60-62˚C for 32 sec. Cotton ubiquitin gene (Gen-Bank EU604080) or Arabidopsis actin gene (GenBank AY064043.1) was used as an endogenous reference for data normalization, and S1 and S3 Tables list the primer sequences used in the present study. The relative expression of the target genes at the mRNA level was determined using the 2 −ΔΔCt method [19]. Relative gene expression was analyzed by one-way ANOVA using the SAS 8.1. The confidence level of all analyses was set at 95%, and P<0.05 was considered statistically significant.

Gene function analysis by VIGS
pCLCrVA and pCLCrVB, the cotton leaf crumple virus (CLCrV)-based vectors, were used for VIGS [11]. RT-PCR was employed to amplify the candidate genes selected from DEGs (primers were listed in S2 Table), and then the amplicons were inserted into pCLCrVA vector. Moreover, pCLCrVA, pCLCrVB and derivatives of pCLCrVA were respectively transferred into Agrobacterium tumefaciens EHA105 via electroporation. The solution of Agrobacterium harboring pCLCrVA or one of its derivatives was mixed with an equal volume of Agrobacterium harboring pCLCrVB. Then, mixture solution was infiltrated into two fully expanded cotyledons of 10-day-old upland cotton seedlings (Zhongzhimian KV3 and 86-1) as previously described [20]. The CLA1 (endogenous cloroplastos alterados 1) gene was used as the control in VIGS, which encoded 1-deoxyxylulose 5-phosphate synthase, the first enzyme of the 2-Cmethyl-D-erythritol-4-phosphate pathway involved in chloroplast development [21]. Each experiment was performed at least in triplicate with more than 20 plants for each construct per repeat. The corresponding primers were used to determine the expression levels of target genes in the leaf tissues of cotton seedling 20 days after VIGS treatment (S3 Table).

Arabidopsis transformation
The coding sequence of GhSKIP35 was cloned into the plant expression vector pPZP111 [22]. The resulting plasmid pPZP111-GhSKIP35 was introduced into A. tumefaciens strain EHA105. Agrobacterium-mediated flower dipping transformation method was used to generate transgenic Arabidopsis plants [23], which were then selected on kanamycin and confirmed by qRT-PCR. Total RNA of wild-type Arabidopsis plants and GhSKIP35 transgenic plants were extracted by TRIzol method. Purified RNA was then reversely transcribed into cDNA. Arabidopsis actin gene (GenBank AY064043) was selected as an endogenous reference gene (primer sequences were listed in S3 Table). qRT-PCR was performed as above mentioned.

Assessment of disease development and statistical analysis
When GhCLA1-silenced plants exhibited the photobleaching phenotype at 20 days after Agrobacterium infiltration, the other target gene silenced plants and negative control plants were subjected to V. dahliae infection by root dipping in conidial suspension (10 7 conidia/mL) for 30 min. In order to assess the disease development, the plants were then replanted in fresh soil. According to the percentage of affected plant tissues, such as chlorosis, necrosis, wilt and defoliation, a Verticillium wilt disease index (DI) scaling from 0 to 4 was introduced to evaluate the disease severity [24]. The DI of Verticillium wilt was determined using the following formula: Data from three separate experiments were analyzed by ANOVA using the SAS 8.1. The confidence level of all analyses was set at 95%, and P < 0.05 was considered as statistically significant.

Variety Zhongzhimian KV3 shows high resistance to V. dahliae inoculation
The highly toxic and defoliant pathogenic V. dahliae strain (V991) was used to inoculate the upland cotton variety Zhongzhimian KV3 (resistant) and 86-1 (susceptible). The disease index (DI) of Zhongzhimian KV3 and 86-1 are 5.303 ± 0.35 and 47.27±2.79 _in the artificial disease nursery, respectively [8,24]. Our data showed significant difference between susceptible and resistant plants in terms of disease symptoms and growth status at 10 days post inoculation (dpi) in the green house. In addition, we reproducibly observed such a process of symptom development after the V. dahliae infection.
Transcriptome characterization of Zhongzhimian KV3 infected with V. dahliae by high-throughput RNA-sequencing To globally identify resistance genes of upland cotton against V. dahliae, we used the highly resistant variety Zhongzhimian KV3 to perform comparative transcriptomics analyses. As the first barrier against V. dahliae invasion, the root tissues inoculated with V. dahliae strain V991 for 24 h and the mock-inoculated root tissues were used to construct cDNA libraries, which were sequenced by Illumina platform. After filtration of raw reads and quality control, about 52 millions of clean reads were generated from each sample, which was sufficient for quantitative analysis of gene expression. To reveal the molecular events behind the gene expression, it is an important step to match the clean reads to genes in order to annotate sequences [25]. In the present study, all clean reads were aligned to the reference genome (http://cgp.genomics. org.cn/page/species/ContactServlet). S4 Table shows the summarized information of the reads.
A total of 31,319 DEGs were detected during the cotton defense response to V. dahliae using RPKM-based method. Moreover, 4,794 DEGs were obtained using FDR 0.01 and the absolute value of Log 2 Ratio!2 as the threshold. Among them, 820 up-regulated genes and 3,974 down-regulated genes were screened, which might be involved in the V. dahliae-G. hirsutum interaction. According to the DEGs expression pattern upon V. dahliae infection, down-regulated model contained the maximum DEGs, showing that negative regulation might play important roles against Verticillium wilt in cotton.

GO clustering and pathway enrichment analysis of DEGs
To better recognize the main biological functions of DEGs, GO was used to comprehensively describe properties of genes, which were divided into various categories, including cellular component (3,025 DEGs), molecular function (2,770 DEGs) and biological process (3,026 DEGs) (Fig 1). 'Primary metabolic process' and 'response to stimulus' were the most highly represented groups in the biological process category, which might be a rapid response of cotton to pathogen. Genes involved in other important biological processes, such as biosynthetic process (10.57%), response to abiotic stimulus (6.8%), signaling (5.53%) and response to biotic stimulus (3.44%), were also identified. In the molecular function category, the two most abundant sub-categories were catalytic activity (35.29%) and binding (31.33%). In the cellular component category, a large number of unigenes were categorized into cell part (29.17%), intracellular part (28.58%) and membrane (13.68%).
In addition, 2,925 DEGs were assigned to 127 KEGG pathways. Table 1

Validation of RNASeq by qRT-PCR
To further validate the reliability of the Illumina sequencing analysis, 18 DEGs were randomly selected for further analysis using qRT-PCR. Their expression trends (Table 2) were almost consistent with the transcriptome sequencing analysis. Among the 18 DEGs, seven DEGs were involved in plant-pathogen interaction pathway and expressed at higher levels than those in non-inoculated plants. Five DEGs were predicted to be involved in plant hormone signal transduction pathway, that jasmonate ZIM domain-containing protein, serine/threonine-protein kinase SAPK1, and protein phosphatase 2 were up-regulated upon V. dahliae infection, while Brassinazole-resistant 1 protein and EIN3-binding F-box protein were down-regulated. Other DEGs, cytochrome P450 71D10 and 3'-N-debenzoyl-2'-deoxytaxol N-benzoyltransferase related to phenylpropanoid biosynthesis pathway, spotted leaf protein, RING finger, CHY zinc finger domain-containing protein 1 and cullin 1 involved in ubiquitin-mediated proteolysis, ribulose bisphosphate carboxylase/oxygenaseactivase 1 involved in proteasome pathway, were all up-regulated in resistant variety Zhongzhimian KV3 compared with susceptible variety 86-1 upon V. dahliae infection (Fig 2).

DEGs involved in plant-pathogen interaction pathway are significantly enriched
In the present study, 295 DEGs were found in 'plant-pathogen interaction' pathway (Ko0426), including 71 up-regulated genes and 224 down-regulated genes. Among them, genes encoding RLKs made up the largest group, which were up-regulated upon V. dahliae infection, including leucine-rich repeat receptor-like serine/threonine-protein kinases (At3g14840, At3g47570), receptor-like protein 12, GbVe1, Verticillium wilt resistance-like protein, receptor protein  MAPKKK (mitogen-activated protein kinase kinase kinase) were up-regulated upon pathogen attack.

WRKYs and RLKs are required for resistance against V. dahliae in cotton
Since VIGS technique has been successfully used to describe gene function in cotton, we employed this approach to further assess the roles of three WRKY genes (WRKY2, WRKY29 and WRKY13) and two RLK genes (LRR receptor-like serine/threonine-protein kinase FLS2, GhFLS2 and G-type lectin S-receptor-like serine/threonine-protein kinase GhGsSRK). First, the cotton CLA1 gene was used as a control to monitor the efficiency of VIGS in cotton. By 20 dpi with CLA1 gene silenced plants, the newly emerging leaves showed extensive chlorosis phenotype, while non-silenced plants (plants infected with CLCrV vector without CLA1 gene) exhibited wild-type phenotype (Fig 3A). These findings suggested that the CLCrV vector could effectively silence the endogenous gene in cotton without obvious effect on plant growth. Next, we respectively cloned five candidate gene fragments into pCLCrVA and infiltrated into two fully expanded cotyledons of 10-day-old cotton seedlings together with pCLCrVB at a ratio of 1:1 by A. tumefaciens, while pCLCrVA empty vector was used as a negative control. At least 20 plants were infiltrated for each construct per repeat. When plants infiltrated with  GhCLA1 displayed the photobleaching phenotype at 20 dpi, the top leaves were harvested to determine the expressions of corresponding genes. Compared with the control, the expressions of all five genes were decreased approximately more than 51.61% (Fig 3B), and then the silenced and negative control plants were subjected to V. dahliae infection. GhFLS2, GhGsSRK, GhWRKY2 and GhWRKY29 were up-regulated upon V. dahliae infection. We respectively silenced these genes in resistant variety Zhongzhimian KV3. GhGsSRK-, GhWRKY2-and GhWRKY29-silenced plants exhibited more severe symptoms compared with the vector control plants (Fig 3C). The average DI of GhGsSRK-, GhWRKY2-and GhWRKY29-silenced plants was 47 Fig 3D). However, GhFLS2-silenced plants did not compromise cotton resistance to V. dahliae (Fig 3C). GhWRKY13 was down-regulated upon V. dahliae infection. We silenced GhWRKY13 gene in susceptible variety 86-1. Data showed that the GhWRKY13 silencing significantly improved the resistance of susceptible variety to V. dahliae compared with the control plants ( Fig 3C). The average DI of GhWRKY13-silenced plants was 29.24±2.13, while that of control (CLCrV-00) and wild-type plants 86-1 was 51.77±1.66 and 42.57±1.43, respectively ( Fig 3D). Taken together, we identified GhGsSRK, GhWRKY2, GhWRKY29 and GhWRKY13 as important components in upland cotton resistance to V. dahliae infection using VIGS assays.

Phenylpropanoid pathway-related genes
In cotton, the phenylpropanoid pathway plays an important role in response to V. dahliae. Cytochrome P450 monooxygenases, including a number of cytochrome P450 (CYP) species, play an important role in the biosynthesis of phenylpropanoids [33]. CYP71D and CYP736 were up-regulated upon V. dahliae infection through RNA-Sequencing. We silenced these two candidate genes in resistant variety Zhongzhimian KV3, respectively (Fig 4C). Severe symptoms were observed in CYP71D-or CYP736-silenced plants compared with vector control plants in response to V. dahliae (Fig 4A). The DI of silenced plants was increased compared with the control plants ( Fig 4B). Shikimate/quinate hydroxycinnamoyl transferase (HCT) functions in the early steps of lignin biosynthesis, and down-regulated HCT in the plants leads to considerable reductions in lignin content [34]. Silenced HCT in the resistant variety Zhongzhimian KV3 resulted in increased susceptibility of resistant cotton to V. dahliae infection (Fig 4A and 4B). The results showed that these three genes were tightly associated with defense response, suggesting that phenylpropanoid synthesis pathway played an important role in cotton resistance to Verticillium wilt. Over-expression of GhSKIP35 contributes to Arabidopsis thaliana resistance to Verticillium wilt It has been reported that ubiquitination plays an important role in leucine-rich repeat (NLR) resistance (R) protein-mediated immunity [35,36]. In cotton, the ubiqutin-protein ligase family is associated with the defense against Verticillium wilt [7]. The SKP1 protein is an evolutionarily conserved subunit of SCF-type E3 ubiquitin ligases that mediate the ubiquitylation of proteins. DEG of SKP1-interacting partner 35 was identified from RNA-sequencing. We amplified the full-length candidate SKP1-interacting partner 35 using RACE-PCR, which was named as GhSKIP35 and contained an open reading frame of 1,863 nucleotides putatively encoding a peptide of 620 amino acid residues [37]. Knockdown of GhSKIP35 resulted in increased susceptibility of resistant cotton to V. dahliae infection (Fig 5B2). An over-expression strategy was used to assess the function of GhSKIP35. Arabidopsis plant was used in the experiment due to technical difficulties and the long duration of cotton transformation. More than 20 transgenic lines constitutively over-expressing GhSKIP35 were obtained. Wild-type and transgenic plants were subjected to V. dahliae infection using root dipping method. Disease symptoms were examined at 15 dpi. Obvious stunting and more necrotic symptoms were found in wild-type plants compared with transgenic plants (Fig 5B1), showing that the spread of necrosis and stunting in transgenic plants was inhibited by GhSKIP35 over-expression. The DI of the transgenic plants was lower than that of the wild-type plants ( Fig 5C). These results suggested that over-expression of GhSKIP35 conferred Arabidopsis plants increased tolerance to Verticillium wilt.

Discussion
Once facing with pathogen attack, plants instigate sophisticated mechanisms in response to pathogen invasion through both constitutive and inducible defenses [38]. In recent years, transcriptome studies have been conducted to reveal the complexity of V. dahliae infection in cotton [1,9,10]. However, the molecular interaction between cotton and V. dahliae remains largely elusive, especially for upland cotton. Currently, the complete gene profiling of cotton species (G. raimondii, G. arboretum and G. hirsutum) is available [39,40,41], and genome database provides convenient for understanding of the cotton defense response to V. dahliae. RNA-Sequencing technology provides huge amounts of first-hand data, which might involve in cotton defense mechanism against V. dahliae infection. Therefore, it is urgent to identify and characterize the functions of these first-hand defense responsive genes on a genome-wide scale, which not only helps understand the resistance mechanism, but also provides support for candidate gene selection to breed disease-resistant varieties by genetic engineering. Because of the difficulty of traditional cotton genetic transformation technique, it is a daunting task to identify the functional roles of DEGs during pathogen attack. With the development of VIGS technology, large-scale functional analysis of individual genes by knocking down the expression of endogenous genes becomes realized.
In the present study, we dissected the transcriptional response of G. hirsutum against V. dahliae infection using the Illumina-based paired-end deep sequencing platform for genomescale analyses. Investigations on the relationships between cotton and V. dahliae revealed that the induction of most defense response occurred early during the interaction [42, 43]. Xu et al [1] selected four time points of V. dahliae inoculation (such as 4, 12, 24 and 48 h) to analyze the defense response by RNA-sequcing, and the results showed that cotton begins to reponse the pathogen infection after 4 hpi. Here, mock-inoculated and 24 h after V. dahliae inoculation were taken for analysis and 52 millions of clean reads were generated from each sample. The data were sufficient for quantitative analysis of gene expression. A total of 4,794 DEGs were either up-or down-regulated by V. dahliae attack. With VIGS technology, we performed functional analysis of several candidate DEGs by suppressing the expressions of endogenous genes in G. hirsutum. The loss-of-function assays provided important evidence for functional genomic studies of cotton genes.
DEGs involved in "plant-pathogen interaction" and 'plant hormone signal' pathways were significantly enriched in the root tissues of resistant G. hirsutum upon V. dahliae infection. Rapid changes in the cytosolic concentration of the secondary messenger Ca 2+ are crucial for downstream responses, such as generation of reactive oxygen species (ROS), activation of calcium-dependent kinase and mitogen-activated protein kinases (MAPKs), and production of defense compounds [27]. A large amount of related DEGs were identified by RNA-sequencing, such as CMLs, CBLs, CaM and EF-hand CBL, which belong to Ca 2+ signals. MAPKK6 and MAPKKK3 belong to MAPK cascades. Respiratory burst oxidase and ferric-chelate reductase were identified and involved in the regulation of ROS status. The SA and jasmonic acid (JA)-ethylene (ET) hormone pathways play a crucial role in the regulation of defence-gene expression [44]. Recent studies also indicated that complex crosstalk among different classes of hormones might modulate the disease resistance, with outcomes dependent on the pathogen lifestyles and the genetic constitution of the host [45]. Eight DEGs genes functions indicated in SA (GhPUB17, GhTGA7 and GhPR1), JA (GhJAZ10 and GhbHLH18), ET (GhEBF1), cytokinine (GhE13L13) and BR (GhBZR1) signal pathways were silenced in upland cotton, and demonstrated that these phytohormones pathway related genes are important components in response to V. dahliae infection [20].
In this study, we identified some important components that are involved in defense response against V. dahliae infection. WRKY gene family is one group, which might play an important role against Verticillium wilt in cotton. Silencing of GhWRKY2 or GhWRKY29 in Zhongzhimian KV3 resulted in a more severe symptom in the silenced plants compared with the vector control plants, indicating that GhWRKY2 and GhWRKY29 largely contributed to cotton resistance to V. dahliae infection. Silencing of GhWRKY13 in susceptible variety 86-1 significantly reduced symptoms compared with the vector control plants, revealing that GhWRKY13 served as a negative regulator of resistance to V. dahliae in upland cotton. As a large gene family, WRKY has received increasing attention due to its roles in plant defense. For example, WRKY29 has been shown to be induced by MAPK cascade in Arabidopsis and tobacco, conferring resistance to both bacterial and fungal pathogens [46,47]. WRKY2 isolated from pepper plays a crucial role in early defense responses to biotic and abiotic stresses [48]. Moreover, the lignin pathway and xylem development can be affected by the grapevine transcription factor WRKY2 in tobacco [49]. In cotton, we also identified two WRKY genes (GhWRKY22, GhWRKY33) from full-length cDNA library, and confirmed that GhWRKY22 and GhWRKY33 were associated with defense response against Verticillium infection [8].
The other group was RLKs genes that were important components involved in "plantpathogen interaction" pathway. Typically, RLKs directly perceive PAMPs by pattern recognition receptors (PRRs) [55], including well characterized LRR-RLKs [56]. Ve1 was important and suggested as being required for resistance against race 1 of V. dahliae [57]. The tomatohomologous Gbve1 gene was up-regulated upon V. dahliae infection. The expression level of Gbve1 was increased 3.27-fold in the inoculated resistant cotton cultivar than that in noninoculated control. It has been found to be involved in defense response against Verticillium [58]. Another important RLKs gene was G-type lectin receptor kinase, GhGsSRK. We silenced GhGsSRK in resistant cotton plants, and the result showed that it played an important role against Verticillium wilt via VIGS analysis. Chen et al. [59] identified a G-type lectin receptor kinase (Pi-d2) from rice, which facilitates resistance against the fungal pathogen Magnaporthe grisea, the causal agent of rice blast.
Other important component was cytochrome P450 (CYPs) gene family. It has been reported that hundreds of CYPs are involved in the phenypropanoid metabolism [60]. Phenylpropanoid pathway was confirmed to play an important role in response to V. dahliae [1,61,62]. CYP71D and CYP736 were up-regulated upon V. dahliae attack. They are involved in ferulate 5-hydroxylase (F5H) synthesis. It has been reported that F5H is cytochrome P450-dependent monooxygenase in phenylpropanoid metabolisms. Over-expression of F5H enhances lignin syringyl contents, which may function in response to various stresses [63,64]. Knockdown of CYP71D and CYP736 in resistant upland cotton variety Zhongzhimian KV3, resulting in the loss of resistance to Verticillium wilt, that might lead to affect F5H synthesis.
In conclusion, our data provided significant insights into the molecular mechanism of cotton responses to V. dahliae, especially a better understanding of early defense response of upland cotton to V. dahliae. We screened several important pathways and components involved in upland cotton defense responses against V. dahliae infection by using VIGS and over-expression stratigies for gene functional characterization. This study provided an important strategy of high-throughput screening and identification of genes involved in defense responses, and our data also offered candidate genes for further genetic engineering of Verticillium wilt tolerance in upland cotton.