Identification of cell wall-associated kinases as important regulators involved in Gossypium hirsutum resistance to Verticillium dahliae

Background Verticillium wilt, caused by the soil borne fungus Verticillium dahliae, is a major threat to cotton production worldwide. An increasing number of findings indicate that WAK genes participate in plant−pathogen interactions, but their roles in cotton resistance to V. dahliae remain largely unclear. Results Here, we carried out a genome-wide analysis of WAK gene family in Gossypium hirsutum that resulted in the identification of 81 putative GhWAKs, which were all predicated to be localized on plasma membrane. In which, GhWAK77 as a representative was further located in tobacco epidermal cells using transient expression of fluorescent fusion proteins. All GhWAKs could be classified into seven groups according to their diverse protein domains, indicating that they might sense different outside signals to trigger intracellular signaling pathways that were response to various environmental stresses. A lot of cis-regulatory elements were predicted in the upstream region of GhWAKs and classified into four main groups including hormones, biotic, abiotic and light. As many as 28 GhWAKs, playing a potential role in the interaction between cotton and V. dahliae, were screened out by RNA-seq and qRT-PCR. To further study the function of GhWAKs in cotton resistance to V. dahliae, VIGS technology was used to silence GhWAKs. At 20 dpi, VIGSed plants exhibited more chlorosis and wilting than the control plants. The disease indices of VIGSed plants were also significantly higher than those of the control. Furthermore, silencing of GhWAKs significantly affected the expression of JA- and SA-related marker genes, increased the spread of V. dahliae in the cotton stems, dramatically compromised V. dahliae-induced accumulation of lignin, H2O2 and NO, but enhanced POD activity. Conclusion Our study presents a comprehensive analysis on cotton WAK gene family for the first time. Expression analysis and VIGS assay provided direct evidences on GhWAKs participation in the cotton resistance to V. dahliae. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-02992-w.


Background
Tetraploid Gossypium hirsutum is the most widely cultivated cotton species in the world and represents an important source of natural fiber and oilseed. Verticillium wilt, caused by the soil borne fungus Verticillium dahliae, is a major threat to cotton production [1]. Identification and characterization of genes associated with resistance is an important basis for potential understanding on the interaction between cotton and V. dahliae, which is necessary for the development of novel disease management methods and new varieties resistant to Verticillium wilt.
Plants live in a complex environment crowded with biotic stresses mainly caused by various phytopathogens and pests, and expose to abiotic stresses including cold, hot, drought and salinity. To overcome these stress challenges, plants have evolved a complex and efficient defense signaling network, which includes monitoring systems to perceive different stress-derived signals triggering specific defense responses [2]. Cell wall, a dynamic structure surrounding plant cell, has emerged as an essential monitoring system [3,4]. Some receptor-like kinases (RLKs) have been identified as cell wall integrity sensors that are responsible for the communication between the cell wall and cytoplasm. Typically, RLKs contain a signal peptide (SP), transmembrane (TM) domain, and cytoplasmic kinase domain. They can be classified into more than 21 subfamilies according to their diverse extracellular domains [5]. Of which, wall-associated kinases (WAKs) are distinguished from the other RLKs by the presence of their unique extracellular epidermal growth factor (EGF)-like domains [5,6].
In Arabidopsis thaliana, WAKs are encoded by 5 WAKs and 22 WAKLs (WAK-like genes) [7]. So far, WAK gene family was also identified in other plants, including Oryza sativa [8], Brassica rapa [9] and Populus trichocarpa [10]. It has been demonstrated that some WAKs are involved in plant development, abiotic and biotic stress responsiveness. Notably, most of WAKs were characterized from Arabidopsis and rice. Arabidopsis AtWAK1, the first identified WAK gene in plant, was shown to contribute to the immune response [11,12]. A rice WAK gene, OsDEES1 (DEFECT IN EARLY EMBRYO SAC1), played a role in the regulation of early embryo sac development [13]. OsiWAK1 (O. sativa indica WAK-1) and HvWAK1 (Hordeum vulgare WAK-1) were involved in plant root development [14,15]. Xa4, encoding a WAK in rice, conferred race-specific durable resistance against Xanthomonas oryzae pv. oryzae by reinforcing the cell wall and increasing the production of jasmonate-isoleucine and phytoalexins [16]. OsWAK1 (O. sativa WAK) and OsWAK25 were upregulated by wounding and salicylic acid (SA), and their overexpression led to higher resistance in transgenic rice lines against Magnaporthe oryzae [17,18]. The other four rice WAK genes, including OsWAK14, OsWAK91, OsWAK92 and OsWAK112d, were also suggested to be required for resistance to M. oryzae by loss-of-function mutants [19]. Beyond rice and Arabidopsis, WAKs have been characterized in response to pathogens as well in other plants, such as tomato SlWAK1 (conferring resistance to Pseudomonas syringae) [20], maize ZmWAK (conferring resistance to Sporisorium reilianum) [21] and ZmWAK-RLK1 (conferring resistance to Setosphaeria turcica) [22].
An increasing number of findings indicate that WAK genes participate in plant−pathogen interactions. Therefore, in our study, we used the latest G. hirsutum genome sequence data (HAU version 1.1 [23] to explore the WAK gene family, representing the first genome-wide identification of GhWAKs. Moreover, two GhWAKs were functionally characterized in response to V. dahliae infection using VIGS (virus induced gene silencing).

GhWAKs identification and localization
In total, 81 GhWAKs as candidates were identified and named according to their chromosomal locations. These GhWAKs were marked on the physical map of 18 chromosomes (Fig. 1a) and one scaffold664 (GhWAK65). A total of 34 and 46 GhWAKs were distributed in the A and D sub-genomes, respectively. Chromosome D02 harbored the largest number of GhWAKs with 20 genes. Six pairs of tandem duplication events were found, including GhWAK16/17, GhWAK36/37, GhWAK43/44-49, GhWAK50/52, GhWAK61/62 and GhWAK69/70/71. These results revealed that the evolution and expansion of GhWAKs happened in G. hirsutum, especially on chromosome D02. The detailed information about GhWAKs, including gene ID, open reading frame (ORF) length, amino acid length, protein molecular weight and isoelectric point, instability index and subcellular localization, was listed in Table 1.
All GhWAKs were predicated to be localized on plasma membrane (PM) ( Table 1). In which, GhWAK77 as a representative was further located in tobacco epidermal cells using transient expression of fluorescent fusion proteins. The images clearly showed that fluorescent signal corresponding to the sole gfp (green fluorescent protein) gene was observed in PM, cytoplasm and nucleus. However, the fluorescent signal corresponding to GhWAK77-gfp was solely shown in PM (Fig.  1b). These suggested that GhWAKs might be a potential connector responsible for communication between inside and outside of the cell.  The majority of GhWAKs have 3-4 introns and show similar exon-intron structure ( Figure S1). A total of six conserved protein domains were identified in GhWAKs, including GUB_WAK_bind (wall-associated receptor kinase galacturonan-binding, PF13947), WAK (wall-associated kinase, PF08488), WAK assoc. (wall-associated receptor kinase C-terminal, PF14380), EGF (EGF, PF00008; cEGF, PF12662; hEGF, PF12661; EGF_CA,  2a). Cytoplasmic, extracellular and TM regions were predicated in the majority of GhWAKs, further indicating that they were PM proteins. Typical WAK encodes a transmembrane protein with a cytoplasmic kinase domain and an extracellular region. However, several proteins showed uncommon structural characteristics, such as the kinase domain in extracellular region, double TMs and kinase domains. All GhWAKs were classified into seven groups according to their protein domain analysis (Fig. 2b)

Prediction of putative cis-regulatory elements in GhWAK promoters
The 2-kb region upstream of the translation start site of all GhWAKs were considered the promoter and analyzed for the potential roles of cis-regulatory elements (Fig. 3). These cis-regulatory elements were classified into four main groups including hormones, biotic, abiotic and light. Twelve hormone-responsive regulatory elements associated with abscisic acid (ABA) (ABRE, ABRE4 and AT-ABRE), auxin (IAA) (AuxRR-core, TGA-box and TGA-element), methyl jasmonate (MeJA) (CGTCAmotif), gibberellin (GA) (GARE-motif, P-box and TATC-box), SA (TCA-element) and ethylene (ET) (ERE), were identified. Of which, ABRE-motif, CGTCAmotif and ERE were enriched in the most of GhWAK promoters, indicating that they might be widely induced by ABA, JA and ET. The biotic stress-related regulatory elements, such as AT-rich, TC-rich repeats, W-box, WUN-motif, WRE3, JERE and box S, were involved in elicitor-mediated activation, wounding and pathogen responsiveness. In addition, eight abiotic-responsive regulatory elements, associated with anaerobic induction (ARE and GC-motif), low-temperature responsiveness (LTR), drought-inducibility (MBS, DRE core and DRE1), heat shock, osmotic stress, low pH, nutrient starvation (STRE) and stress-related (TCA), were identified in the GhWAK promoter regions. Moreover, various light-responsive elements were present in the promoters of GhWAKs. Especially, Box 4 and G-Box were widely harbored. These results indicated that GhWAKs might play vital roles in the response to various stresses, hormones and light.

GhWAKs were significantly induced by V. dahliae infection
To identify GhWAKs that were related to V. dahliae infection, two-fold changes were applied in transcript expression profiles from RNA-seq as minimum cutoffs. As a result, 26 GhWAKs were screened out, including 17 upregulated and 9 down-regulated genes (Fig. 4a). Of which, 11 GhWAKs, including GhWAK5, GhWAK9, GhWAK77, GhWAK10, GhWAK45, GhWAK47, GhWAK78, GhWAK48, GhWAK31, GhWAK26 and GhWAK72, were significantly up-regulated in at least three time points, suggesting that they continuously responded to V. dahliae infection. Their expression profiles were further verified through real-time quantitative reverse transcription PCR (qRT-PCR). The expression results of them in response to the V. dahliae infection from qRT-PCR were consistent with those found in RNA-seq data (Fig. 4b). Due to the high degree of sequence similarity in GhWAKs family, it was difficult to design specific primers for four gene pairs, including GhWAK4/GhWAK45, GhWAK5/GhWAK49, GhWAK10/GhWAK55, and GhWAK31/GhWAK77. The results of qRT-PCR indicated that these four pairs of GhWAKs were dramatically up-regulated. According to RNA-seq data, GhWAK4, GhWAK49 and GhWAK55 did not show to be up-regulated. Thus, the expression changes found using qRT-PCR probably more represent the responses of GhWAK45, GhWAK5 and GhWAK10 to V. dahliae infection. In addition, the other 45 GhWAKs that did not show differential expression in RNA-seq data were further detected through qRT-PCR. As a result, GhWAK1 and GhWAK69 showing higher transcription levels in cotton seedlings inoculated with V. dahliae than that in control was screened out complementally (Fig. 4c). Finally, a total of 28 GhWAKs were found to play a potential role in the interaction between cotton and V. dahliae.
Silencing GhWAKs compromised cotton resistance to Verticillium wilt GhWAK26 and GhWAK77 showed obviously and persistently up-regulated expression to the infection from Fig. 3 Potential cis-elements in a 2 kb 5′ flanking region upstream from the start codon of each GhWAK. The number of each cis-element was shown, and the back-color changes from blue to red as the number increase. All cis-regulatory elements were classified into four groups, including hormones, biotic, abiotic and light V. dahliae (Fig. 4). In addition, they contain cis-elements in their promoters associated with MeJA and SA, which play key roles in cotton resistance to V. dahliae. Thus, to further reveal the function of GhWAKs in cotton resistance to V. dahliae, GhWAK26 and GhWAK77 were prioritized for study as representatives here using tobacco rattle virus (TRV) based VIGS system. At approximately two weeks post-infiltration with a mixture of Agrobacterium cultures containing pTRV1 and pTRV2-CLA1, a strong photobleaching phenotype was shown on . FC is a ratio of treatment FPKM to control FPKM that obtained from the RNA-Seq data. Higher and lower transcriptional level are indicated by pink and blue, respectively, and no detected expression is indicated by dark color. b and c, The expression analysis of GhWAKs from G. hirsutum ND601 induced with V. dahliae by qRT-PCR. The relative gene expression level was calculated using the comparative 2 -ΔΔCt method with GhHis3 as internal control. The bar represents the standard error (SE) calculated from three independent experiments. Asterisks indicate statistically significant differences according to Sidak's multiple comparisons test (**P < 0.001, ***P < 0.01; ns, no significant) the newly emerging true leaves (Fig. 5a), indicating that VIGS system worked well. Then, the expression of GhWAK26 and GhWAK77 was detected in the leaves infiltrated with pTRV2-GhWAK26 and pTRV2-GhWAK77, respectively. As shown in Fig. 5b, the expression of GhWAK26 and GhWAK77 was reduced by about 80%, suggesting VIGS triggered their silencing in cotton plants. At 20 days post inoculation (dpi), VIGSed plants ( Fig. 5d and e) exhibited more chlorosis and wilting than the control plants infiltrated with Agrobacterium cultures containing empty vector pTRV1 and pTRV2 (Fig. 5c). The disease indices of VIGSed plants were also significantly higher than those of the control at 15 dpi and 20 dpi (Fig.  5f). Therefore, the results of VIGS assays suggested that GhWAK26 and GhWAK77 were important participants in cotton resistance to V. dahliae infection.
Silencing GhWAKs increased the spread of V. dahliae in cotton stems After inoculation, V. dahliae in cotton stem was detected by PCR. No specific amplification products from V. dahliae were shown in CK at 5 dpi and 7 dpi, indicating that V. dahliae had not yet invaded the stems or multiplied in large quantities ( Fig. 6a and Figure S2). However, at 5 dpi, few specific products from V. dahliae were amplified in GhWAK26-silenced and GhWAK77-silenced plants, representing a small amount of pathogen invasion. Further, at 7 dpi, the bright bands amplified from GhWAK26-silenced and GhWAK77-silenced plant stems appeared on agarose gels, indicating that V. dahliae had invaded largely. In addition, pathogen isolation on potato dextrose agar (PDA) showed that a large number of V. dahliae grew out from the stems of GhWAK26-silenced and GhWAK77-silenced cotton plants, while no mycelium was shown from the control (Fig. 6b). Both PCR detection and PDA culture results suggested that silencing GhWAKs significantly increased the spread of V. dahliae in the cotton stems.
Lignin is considered to play an important role in preventing cotton from the infection of V. dahliae. Therefore, we further compared the changes of lignin content in GhWAK-silenced cotton stems with CK. The results showed that the lignin content in GhWAK-silenced plants was significantly lower than that in CK (Fig. 6c), which might affect the stem structure and then reduce the prevention of cotton from V. dahliae infection. The content of H 2 O 2 and NO, and POD activity in GhWAK-silenced plants inoculated with V. dahliae were further measured. GhWAKs silencing caused lower levels of H 2 O 2 at 6 h post inoculation (hpi), 12 hpi and 24 hpi ( Fig. 7a and b). Both GhWAK26and GhWAK77-silenced plants accumulated greatly depressed levels of NO comparing with CK ( Fig. 7c and d). However, the activity of POD significantly elevated in GhWAK26and GhWAK77-silenced plants at 6 hpi, 24 hpi and 48 hpi, except at 12 hpi ( Fig. 7e and f).
Silencing GhWAKs significantly affected the expression of JA and SA-related marker genes Further, the expression of several JA and SA-related marker genes involved in plant defense signaling pathways was detected. The expression of JAZ1 (jasmonatezim-domain protein), JAZ3, JAZ6, LOX1 (lipoxygenase) (JA-related marker genes), PR3 (pathogenesis related protein) and NPR1 (nonexpresser of PR protein) (SA-related marker genes) were significantly down-regulated after silencing GhWAK26 in cotton (Fig. 8a). In GhWAK77-silenced plants, JAZ6 and three important genes involved in the SA signaling pathway, including ICS1 (isochorismate synthase), NPR1 and EDS1 (enhanced disease susceptibility), were down-regulated comparing with control. On the contrary, the expression of JAZ1 and LOX1 were significantly up-regulated due to the silencing of GhWAK77 (Fig. 8b). These results indicated that GhWAK26 and GhWAK77 might involve in cotton resistance to V. dahliae through SA and JA signaling pathways.
At the PM, RLKs as cell-surface receptors can perceive and process extracellular danger signals to trigger plant defense responses [28]. WAK belongs to RLK subfamily. All GhWAKs contain a typical eukaryotic kinase domain that is mostly present in intracellular region and relatively well conserved (Fig. 2a). In addition, GhWAKs locate on PMs in all probability (Table 1, Fig. 1), suggesting that GhWAKs have potential roles in communicating between inside and outside of the cell. In order to penetrate plant roots to gain access to the xylem and to spread in the vascular system, V. dahliae usually secretes various toxins and carbohydrate active enzymes, including glycoproteins and cell wall-degrading enzymes [29,30]. Therefore, it is conceivable that V. dahliae infection affects plant cell wall integrity (CWI) and generates some degradation products, which are important defense signals [31]. In the extracellular region, GhWAKs contain five different domains (Fig. 2a), which may sense CWI or interact with different components of these extracellular matrix, such as glycine-rich protein, pectin and oligogalacturonides (OGs) [32][33][34].
At present, the molecular mechanism of WAK-mediated resistance remains largely unknown. However, some defence responses associated with WAKs have been reported, including cell wall reinforcement [16], pathogenesis-related genes activation [18], SA or JA accumulation [27], POD and superoxide dismutase activities [27], and reactive oxygen species (ROS) homeostasis [27]. Here, silencing GhWAKs resulted in Fig. 7 Silencing of GhWAKs dramatically compromised V. dahliae-induced accumulation of H 2 O 2 (a and b) and NO (c and d), but enhanced POD activity (e and f). The results from three biological replicates are shown with mean ± SE. Asterisks represent P values (**P < 0.001, ***P < 0.001; ns, no significant; Sidak's multiple comparisons test) the up-or down-regulation of several genes (Fig. 8) and depressed cotton resistance to V. dahliae. Among them, JAZ and LOX are associated with JA-mediated defense responses [35]. NPR1, ICS1 and EDS1 are associated with SA-mediated defense responses [36]. The two phytohormones, JA and SA, have been known to be involved into the regulation of plant resistance against V. dahliae [37,38]. In addition, some hormone-responsive and biotic stress-related regulatory elements were enriched in the promoters of GhWAKs (Fig. 3). Thus, these findings suggest that GhWAK function as a mediator to active intracellular SA and JA signaling pathways to regulate cotton resistance.
V. dahliae is a vascular pathogen that penetrates the host roots and then extends to other overground parts of plant through the process of transpiration [29,37]. The improvement of physical, chemical and structural barriers, such as ROS, NO, cell wall, lignin, callose and dahliae. JAZs, LOX1 and PR3 are the marker genes involved in JA signaling pathway. ICS1, NPR1 and EDS1 are the marker genes involved in SA signaling pathway. The results from three biological replicates are shown with mean ± SE. Asterisks represent P values (*P < 0.05; **P < 0.01; ***P < 0.001; ns, no significant; Sidak's multiple comparisons test) POD, contributes to preventing expansion and reducing colonization of V. dahliae in cotton tissues [37,[39][40][41]. In this study, more V. dahliae was detected in GhWAK26-silenced or GhWAK77-silenced plants with lower lignin contents than in CK (Fig. 6). Moreover, silencing of GhWAKs in cotton plants dramatically compromised V. dahliae-induced accumulation of H 2 O 2 and NO, but enhanced POD activity (Fig. 7). These findings demonstrate that GhWAKs play roles in preventing pathogen spreading at least in part by regulating the accumulation of lignin, H 2 O 2 and NO, and the activity of POD. Overall, these results augment our knowledge about cotton WAK gene family, and particularly promote the understanding on their function in disease resistance.

Conclusions
In this study, we carried out a genome-wide analysis of WAK gene family in G. hirsutum with the identification of 81 putative GhWAKs, which might sense different outside signals to trigger intracellular signaling pathways that response to various environment-stresses. Of which, 28 GhWAKs with potential roles in the interaction between cotton and V. dahliae were screened out. Silencing GhWAKs could significantly affect the expression of JA-and SA-related marker genes, increased the spread of V. dahliae in the cotton stems, dramatically compromised V. dahliae-induced accumulation of lignin, H 2 O 2 and NO, but enhanced POD activity. These results provided direct evidences that GhWAKs participate in the cotton resistance to V. dahliae. Finally, a model for how GhWAKs were involved in cotton resistance to V. dahliae was proposed (Fig. 9).

Analysis of chromosomal location, genes structure and cis-elements
The information about physical chromosomal locations and gene structures of GhWAKs was extracted from the gene annotations in gene feature format (GFF) files, which were downloaded from the CottonFGD website and analyzed by TBtools software [42]. The potential promoter sequences, 2 kb upstream of GhWAKs, were also extracted from G. hirsutum genome database. The cis-elements in the potential promoters were predicted using PlantCARE databases (http://bioinformatics.psb. ugent.be/webtools/plantcare/html/).

Plant materials and V. dahliae inoculation
The seeds of Nicotiana benthamiana and G. hirsutum cv. Nongda 601 (ND601) were preserved at the State Key Laboratory of North China Crop Improvement and Regulation, Hebei Agricultural University, China. N. benthamiana was grown in the greenhouse about 5 weeks at 21°C with 14/10 h (light/dark) photoperiod. ND601 were grown in the greenhouse at 25°C under a 14-h light/10-h dark cycle with relative humidity about 70%. Cotton seedlings inoculation with V. dahliae strain Linxi 2-1 (10 7 spores ml − 1 ) was performed as previously described [39].

Proteins subcellular localization
The ORF of GhWAK77 (without the stop codon) was amplified by PCR with primers gWAK77-F and gWAK77-R (Table S1), and then introduced into entry vector pDONR™207 by attB/attP recombination reaction, as described by the manufacturer (Invitrogen). The GhWAK77 fragment was transferred from the entry clone to expression vector pEarlyGate103 [43] with attL/ attR recombinant reaction, as described by the manufacturer (Invitrogen). The recombinant expression vector was introduced into Agrobacterium tumefaciens GV3101, cultured and infiltrated into four-week-old tobacco leaves via the method described by [44]. After 2 days, GFP signal in the tobacco leaf epidermal cells was examined using a laser scanning microscope (FluoView FV1000; Olympus).

RNA-seq data and qRT-PCR analysis
The transcription patterns of GhWAKs in cotton roots after inoculation with V. dahliae were analyzed using high-through RNA-seq data published previously [37]. Log 2 Fold change were calculated from FPKM (fragments per kilobase of exon model per million mapped) and used for the heat map of hierarchical clustering with the TBtools v0.67 software [42]. Total RNA was extracted using EASYspin Plant RNA kit (Aidlab, Beijing, China) according to the manufacturer's instructions. The quality and concentration of RNA were detected by 1.5% agarose gel electrophoresis and NanoDrop™ 1000 spectrophotometer (Thermo Fisher Scientific), respectively. cDNA was synthesized with a reverse transcription kit (ReverTra Ace® qPCR RT Master Mix with gDNA Remover, TaKaRa, Dalian, China). qRT-PCR was performed using 7500 Real Time PCR System (Applied Biosystems, USA) with THUNDERBIRD®SYBR® qPCR Mix (TaKaRa, Dalian, China). The 2 -ΔΔCt method was used to calculate the relative expression of genes. GhHis 3 was used as internal reference. Three biological repeats were taken for each treatment.

VIGS assays in cotton
The vectors for VIGS, pTRV1 and pTRV2, were kindly provided by Professor Liu Yule of Tsinghua University [45]. The fragments from GhWAKs were amplificated by PCR and inserted into the pTRV2 vector between EcoR I and Kpn I. The constructed vectors were separately transferred into A. tumefaciens strain GV3101 by freezethaw method [46]. VIGS in cotton was performed as described previously [47]. At least 30 plants were used per treatment, and each treatment was repeated three times. Plant resistance to V. dahliae was assayed by analyzing disease index [48].
Detection and isolation of V. dahliae in cotton stems At 5 dpi and 7 dpi, 1 cm and 0.5 cm of samples excised at a height of 0.5 cm stem above ground were used for detection and isolation of V. dahliae, respectively. V. dahliae detection by PCR was performed using primers P1 and P2 [49]. V. dahliae isolation from cotton stems was carried out according to the previous method [50]. Twenty-four individual plants were sampled for each treatment and repeated three times.

Measurements of NO, H 2 O 2 and POD activity
The first true leaves of cotton seedlings were powdered in the mortar with liquid nitrogen and homogenized using 50 mM sodium phosphate buffer (pH 7.0). After centrifugation (14,000 g, 20 min), the supernatants were used for the determination of NO, H 2 O 2 and POD activity with commercialized assay kits (Nanjing Jiancheng Bioengineering Institute, China), following the manuals. The total protein concentration of the supernatants was measured using Pierce™ BCA Protein Assay Kit (Thermo Scientific).

Primers and statistical analysis
All primers used in this study were listed in Table S1. Differences between measured values were analyzed using software GraphPad Prism® 8 (GraphPad, San Diego, CA, USA). A two-way ANOVA with multiple comparisons (Sidak's test) was used to compare gene expression in cotton roots between inoculated with V. dahliae and inoculated with water (CK) at the same hpi, disease indices, H 2 O 2 and NO content, and POD activity between GhWAK-silenced plants and CK. A one-way ANOVA with Dennett's multiple-comparisons test was used to compare lignin content between GhWAK-silenced plants and CK. The P-value less than 0.05 was assumed to be statistically significant.