Phosphoproteomics Profiling of Cotton (Gossypium hirsutum L.) Roots in Response to Verticillium dahliae Inoculation.

Verticillium wilt is a plant vascular disease causing severe yield and quality losses in many crops and is caused by the soil-borne plant pathogenic fungus Verticillium dahliae. To investigate the molecular mechanisms of the cotton-V. dahliae interaction, a time-course phosphoproteomic analysis of roots of susceptible and resistant cotton lines in response to V. dahliae was performed. In total, 1716 unique phosphoproteins were identified in the susceptible (S) and resistant (R) cotton lines. Of these, 359 phosphoproteins were significantly different in R1 (1 day after V. dahliae inoculation) vs R0 (mock) group and 287 phosphoproteins in R2 (3 days after V. dahliae inoculation) vs R0 group. Moreover, 263 proteins of V. dahliae-regulated phosphoproteins were significantly changed in S1 (1 day after V. dahliae inoculation) vs S0 (mock) group and 197 proteins in S2 (3 days after V. dahliae inoculation) vs S0 group. Thirty phosphoproteins were significantly changed and common to the resistant and susceptible cotton lines following inoculation with V. dahliae. Specifically, 92 phosphoproteins were shared in both in R1 vs R0 and R2 vs R0 but not in susceptible cotton lines. There were 38 common phosphoproteins shared in both S1 vs S0 and S2 vs S0 but not in resistant cotton lines. GO terms and Kyoto Encyclopedia of Genes and Genomes pathway analyses displayed an abundance of known and novel phosphoproteins related to plant-pathogen interactions, signal transduction, and metabolic processes, which were correlated with resistance against fungal infection. These data provide new perspectives and inspiration for understanding molecular defense mechanisms of cotton roots against V. dahliae infection.

Introduction

Cotton (Gossypium spp.) is one of the most important economic crops in the world and is the main source of natural textile fibers. The plant pathogenic fungus Verticillium dahliae Kleb. causes plant vascular diseases[1] and often results in severe cotton yield reductions and losses.[2,3] The pathogen can adhere to the root surface, and its hyphae penetrate the roots and colonize the cortex. They then extend into the xylem and form conidia.[4] Meanwhile, the toxic, cell-wall-degrading enzymes and elicitor-like substances are produced by V. dahliae to suppress plant innate immune responses.[3] It has been found that the tolerant and resistant plant genotypes are less sensitive to V. dahliae and activate defense responses to V. dahliae.(5−7) Although a number of studies report on plant–V. dahliae interactions, there remains limited knowledge about the molecular mechanism of plants against V. dahliae.(5,8,9)

Post-translational modifications play a key role in many cellular processes and regulate protein activity, targeting, or interactions with other proteins.[10] The reversible phosphorylation of proteins plays an important role in the regulation of almost all life activities since it plays a vital role in intracellular signal transduction[11,12] and is involved in regulating cell cycle progression, development, response, differentiation, transformation, and metabolism.[13] Protein kinases are involved in the reversible process (phosphorylation/dephosphorylation), including the transfer and removal of the phosphate groups, which are important signaling mechanisms in eukaryotic cells.[14−18] These mechanisms have been confirmed to play important roles in response to a variety of abiotic and biotic stresses, including cold stress, drought stress, salt stress, and pathogen infections. Deciphering molecular mechanisms of defense in response to pathogen stresses will improve our understanding of the processes in plants.[19,20]

In recent years, large-scale protein phosphorylation analyses have been performed on many plant species with technological advancements in peptide enrichment, high-resolution mass spectrometry (MS), and associated bioinformatics.[20−24] In maize, a total of 4361 phosphorylated protein sites were quantified, of which were significantly regulated upon drought stress. In addition, it is as major functional groups that signaling and transport proteins were identified in the post-translation modifications. Similarly, most of the identified proteins in maize roots grown under high and low iron conditions were classified as signaling proteins. Of these, BAK1, PTI1, AIR12, and 14-3-3 proteins were significantly changed by Fe stress, and PTI1 levels were increased in response to high and low iron stresses.[25] Protein phosphorylation and dephosphorylation are involved in multiple signaling pathways. In plants, phosphorylation and dephosphorylation also play a pivotal role in ABA signaling pathway and mediate the regulation of stomatal movement in response to ABA.[26] The phosphorylation of CIPK is an important factor in Ca2+-mediated signaling, which plays a key role in signal transduction.[27]

To investigate the molecular mechanisms of the cotton–V. dahliae interaction, we analyzed changes of phosphoprotein profiles in the roots of susceptible and resistant cotton genotypes inoculated with V. dahliae using liquid chromatography-tandem mass spectrometry (LC-MS/MS). We also performed proteome-wide label-free quantification using MaxQuant software with intensity-based absolute quantification (iBAQ). This study provides valuable data and novel insights for further investigations into molecular events involved in cotton defense responses to V. dahliae.

Results and Discussion

Identification and Differential Analysis of Cotton Roots in Response to V. dahliae Infection

The strategy of phosphoproteomic profiling from cotton roots inoculated with V. dahliae is shown in Figure S1. A total of 5405 unique phosphosites and 2522 unique phosphopeptides from 1716 unique proteins were identified from the roots of NIL-R and NIL-S (Table S1). Of the 5405 phosphorylation sites, 4702 (87.00%) were from phosphoserine, 619 (11.45%) were from phosphothreonine, and 84 (1.55%) were from phosphotyrosine. Of the 2522 unique phosphopeptides, the singly phosphorylated were very few, only 1.07%, whereas most were doubly or multiply phosphorylated (46.75 and 52.18%, respectively) (Table ).

Table 1

Number of Phosphorylation Sites, Phosphopeptides, and Phosphoproteins Identified in This Study

listing	number
phosphorylation sites	5358
phosphorylation sites (tyr/hr/ser)	84:614:4660
phosphopeptides	2522
phosphopeptides (single/double/multiple)	27:1179:1316
phosphoproteins	1716

Among those 1716 phosphoproteins, differences between the mock and inoculated samples were assessed with a t-test. A threshold p value of <0.05 and the absolute fold change value of phosphorylation of >2 were used as criteria to identify phosphoproteins with a differential abundance. For a total of 639 phosphoproteins, differences between the mock and inoculated samples in cotton roots were assessed. Of these, there were 359 proteins with significant phosphorylation at 1 dpi (R1) compared with the mock (R0) for the resistant line NIL-R [258 proteins were increased and 101 proteins were decreased (Table S2 and Figure A)]. Three days after inoculations (R2), there were 287 proteins with significant phosphoprotein changes compared to the mock (R0) from NIL-R [209 proteins were increased and 78 proteins were decreased (Table S2 and Figure B)]. In the susceptible line NIL-S, there were 263 proteins with significant phosphorylation changes at 1 dpi (S1) compared to the mock (S0) from NIL-S; 144 proteins were increased and 119 proteins were decreased (Table S2 and Figure C). Three days after inoculation (S2), 197 unique proteins had significant changes in phosphorylation between the inoculated and the mock (S0) samples from NIL-S; 121 proteins were increased, whereas 76 proteins were decreased (Table S2 and Figure D).

Figure 1

Hierarchical clustering representation of the differentially expressed phosphorylated peptides from cotton (Gossypium hirsutum) roots in NIL-R (R) and NIL-S (S) lines inoculated with V. dahliae. R0 and S0: mock; R1 and S1: after 1 day of V. dahliae inoculation; R2 and S2: after 3 days of V. dahliae inoculation. (A) R1 vs R0 comparison; (B) R2 vs R0 comparison; (C) S1 vs S0 comparison; and (D) S2 vs S0 comparison. X axis: the hierarchical clusters for samples. Y axis: the hierarchical clusters for phosphorylated peptides. Hierarchical clustering analysis was performed on the log 2 expression of differentially expressed modified peptides displaying similarity among treatments. Red represents significantly increased modified peptides. Green represents significantly decreased modified peptides. Gray represents the modified peptides having no quantitative information.

Phosphopeptide Motif Identification in Cotton Roots in Response to V. dahliae

Protein post-translational modification was performed by an upstream enzyme to identify substrate-specific motifs. Sequence conservation at phosphosites was evaluated for protein post-translational modification and phosphorylation site prediction. Motif-x was used to predict motifs by Motif-X algorithm[28] in this study. In total, 37 Ser motifs (Table S3, sheet 1) and 4 Thr motifs (Table S3, sheet 2) were significantly enriched. Based on the published information on the association of these motifs with protein kinases,[29] 21 novel phosphorylation motifs were identified, including 19 pSer motifs and 2 pThr motifs.

For phosphorylation motifs identified, we found eight corresponding motifs in the PhosphoMotif Finder database (Table S3, sheet 3). Two motifs, .......S..E.... and .......SD.E...., are potential substrates for casein kinase II (CKII), which takes part in the cell-cycle regulation, rhythm regulation, DNA repair, and other physiological regulatory pathways.[30] Three motifs, .......SP......, .......SSP....., and.......TP, belong to proline-directed motifs, which are potential substrates for pro-directed kinase, cyclin-dependent kinase (CDK5), mitogen-activated protein kinase (MAPK), glycogen synthase kinase 3 (GSK-3), and extracellular regulated protein kinases (ERK1/2). Some kinases, such as protein kinase A (PKA), Ca2+/calmodulin-dependent protein kinase (CaM), and calcium-dependent protein kinase (CDPK) could recognize the motifs ....R..S........ These kinases may play roles in signal recognition and transduction of cotton in response to infection by V. dahliae.

GO Annotation and Functional Classification of Significantly Changed Phosphoproteins

To further understand the functional distribution of the significant differentially expressed phosphoproteins identified from cotton roots inoculated with V. dahliae, the phosphoproteins were annotated and classified into cellular component (CC), molecular function (MF), and biological process (BP) categories by Blast2GO program based on level 2 (Figure ). CC analysis showed that protein phosphorylation events occurred in various types of cells, organelles, membranes, macromolecular complexes, and membrane-enclosed lumen both in the resistant and susceptible genotypes (Figure and Table S4). These results are consistent with previous findings that the phosphoproteins were more localized in intracellular organelles, cytoplasm, ribonucleoprotein complexes, etc.[29]

Figure 2

GO annotation of the proteins containing differentially expressed phosphorylated peptides in NIL-R (R) and NIL-S (S) roots inoculated with V. dahliae. X axis: GO terms in three categories: biological processes (marked with red), molecular functions (marked with purple), and cellular compartments (marked with orange); Y axis on the left: the percentage of annotated proteins in all of the differentially expressed proteins; Y axis on the right: the number of the annotated differentially expressed proteins.

BP category provides information on the phosphoproteins over-represented in metabolic processes, biogenesis, biological regulation, or response to biotic/abiotic stress, cellular processes, cellular component organization, and single-organism processes (Figure and Table S4). Furthermore, level 3 scans indicated that 9, 6, 7, and 6 protein kinases were linked to the cotton response to fungal infection in the R1 vs R0, R2 vs R0, S1 vs S0, and S2 vs S0 comparisons, respectively (Table S4). These phosphoproteins are involved in plant hormone signal transduction (EIN2),[31] respiratory burst oxidase (RBOH),[32] plant–pathogen interaction (RIN4),[33] and translation initiation factor (EIF4G)[34] and may be associated with the cotton against V. dahliae.

MF category reveals the over-representation of phosphoproteins involved in catalytic activity and binding. The latter includes enzyme regulator activity, transporter activity, structural molecule activity, and nucleic-acid-binding transcription factor activity (Figure ). In the R1 vs R0 comparison, there were 208 binding proteins and 108 catalytic proteins; the corresponding values were 163 and 74 in the R2 vs R0 comparison, 138 and 83 in S1 vs S0 comparison, and 116 and 55 in S2 vs S0 comparison (Figure S2 and Table S4). A level 5 scan identified 33, 16, 27, and 15 protein kinases in the R1 vs R0, R2 vs R0, S1 vs S0, and S2 vs S0 comparisons, respectively (Table S4). Most of these kinases are essential signaling proteins and receptors for and the tricarboxylic acid (TCA) cycle in plants, such as mitogen-activated protein kinase (MAPK1_3),[35] calcium-dependent protein kinases (CPK),[36] phosphoenolpyruvate carboxykinase (pckA),[37] and protein kinase (PRPF4B, CKI, SRPK3, STK23, PRKWNK, and CSNK1E). These proteins and relevant signaling pathways may play crucial roles in the cotton–V. dahliae interaction.

Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway and Regulatory Roles of Phosphoproteins with Significant Changes

KEGG online service tool KAAS maps KEGG pathways of those proteins with significant differences in phosphorylation. These phosphoproteins are found to be involved in many essential biological pathways, including plant–pathogen interaction, plant hormone signal transduction, HIF-1 signaling pathway, ABC transporters, mTOR signaling pathway, MAPK signaling pathway, PI3K-Akt signaling pathway, etc. (Figure S3 and Table S5). Further research is needed to understand the precise role of these pathways in the cotton against V. dahliae.

Four phosphoproteins related to plant–pathogen interactions were identified: calcium-dependent protein kinase (CDPK28, Gh_A11G1615), respiratory burst oxidase (RBOH, Gh_D05G2471), calcium-binding protein CML (CML41, Gh_A11G1729), and RPM1-interacting protein 4 (RIN4) (Gh_D08G2088) (Figure and Table S5). Compared to that of the mock, the phosphorylation level of CDPK increased significantly at 1 dpi, with 5.5-fold and 4-fold changes observed in the susceptible and resistant cultivars, respectively, but these levels decreased at 3 dpi. The transcript levels for CDPK were 3-fold higher in the roots inoculated with Vd080 than in the mock, 12 hpi, and then decreased thereafter. In contrast, the transcript levels for RBOH were lower in the inoculated roots than in the mock. Both CML41 and RIN4 were not phosphorylated in the cotton roots before inoculation but were significantly phosphorylated after inoculation. The phosphorylation level of RIN4 was 6-fold and 2.7-fold higher in roots of the resistant cultivar than in the resistant cultivar, 1 and 3 dpi, respectively. The transcript levels for CML41 and RIN4 were higher (from 12 hpi onward) in the inoculated roots than in the mock (Figure ). RIN4 is related to plant–pathogen interactions. For example, RIN4 phosphorylation was induced by AvrB (Pseudomonas syringae) targeting a host guardee protein RIN4 in Arabidopsis. Then, phosphorylated RIN4 activates the immune receptor RPM1 to mount defense.[38−41]

Figure 3

Phosphoproteins identified in NIL-R (R) and NIL-S (S) roots inoculated with V. dahliae related to plant–pathogen interaction pathway in the KEGG database.

Figure 4

(A) Protein phosphorylation levels and (B) gene expression analyses of GhCDPK28, GhRBOH, GhCML41, and GhRIN4 in cotton roots.

Many phosphoproteins involved in plant defense signaling pathways were identified in cotton roots. The mTOR signaling pathway, the HIF-1 signaling pathway, and the PI3K-Akt signaling pathway have been usually described in animals and microbes.[42−44] These are related to glycogen metabolism and cell-cycle progression, including reactive oxygen species generation. The small subunit ribosomal protein S6e (RPS6) was identified in the three signaling pathways in this study (Figure S4A). The RPS6 belongs to TIR-NBS-LRR protein, which regulates the resistance of Arabidopsis to P. syringae by effector HopA1.[45] In barley, resistance locus RPS6 as an intermediate enhanced host resistance to wheat stripe rust.[46] These results showed that phosphorylation proteins and the associated signaling pathways play a regulatory role in the cotton–V. dahliae interaction.

Two phosphoproteins were identified in the MAPK signaling pathway: MAP kinase substrate 1 (MKS1) and mitogen-activated protein kinase 1/3 (MAPK1/3). Also, this is indicated that the phosphorylation of the two proteins may be involved in cell death, defense response, and camalexin synthesis. Phosphorylation events were also found in the ethylene (ET) signaling pathway, including the ethylene receptor EIN2 (Figure S4B), which is a positive regulator of ethylene signaling.[47] Another phosphoprotein from the ARR-B family (two-component response regulator) involved in camalexin synthesis (Figure S4C) and associated with cell division and shoot initiation was identified.

Comparison of Phosphoproteins between Resistant and Susceptible Lines

Thirty phosphoproteins were significantly changed common to all of the following comparisons: R1 vs R0, R2 vs R0, S1 vs S0, and S2 vs S0 (Figure and Table S6). Two phosphoproteins (Gh_A05G3831 and Gh_D07G0502) belonging to the adenosine 5′-triphosphate (ATP)-binding cassette were increased in both susceptible and resistant genotypes after inoculation. Another phosphoprotein (Gh_A10G0583) belonging to the ABC transporter G family member 36 was also increased in both lines after inoculation (Table S2). One calcium ion binding protein (Gh_A01G0834) was increased in the resistant lines but was decreased in the susceptible lines. When plants sense and respond to stimuli, including biotic and abiotic stresses, a series of PRRs (pattern recognition receptors) localized in the plasma membrane recognize the damage-associated molecular patterns (DAMPs) and microbe-associated molecular patterns (MAMPs), and induce a rapid [Ca2+]cyt elevation in plant cells.[48,49] The burst of reactive oxygen species and synthesis of defense-related plant hormones, such as jasmonic acid (JA), ethylene (ET), and salicylic acid (SA), are regulated by a series of Ca2+ signaling and proteins relays.[50−52]

Figure 5

Venn diagram representing overlapping of significantly changed phosphoproteins shared between phosphoproteomes of R1 vs R0, R2 vs R0, S1 vs S0, and S2 vs S0 comparisons. R: NIL-R lines; S: NIL-S lines; R0 and S0: mock; R1 and S1: after 1 day of V. dahliae inoculation; R2 and S2: after 3 days of V. dahliae inoculation.

There were 92 phosphoproteins common to both the R1 vs R0 and R2 vs R0 comparisons but not to the S1 vs S0 and S2 vs S0 comparisons. However, there were 38 phosphoproteins common to both the S1 vs S0 and the S2 vs S0 comparisons but not to the R1 vs R0 and R2 vs R0 comparisons (Figure and Table S6). Some of these phosphoproteins were shown to play crucial roles in plant–pathogen interactions. For example, the high mobility group box domain (Gh_A02G0873) was increased (by 5-fold) in resistant lines but not in susceptible lines (Table S2). The HMGB3 protein is a novel plant DAMP, which induces hallmark innate immune responses through receptor-like kinases BKK1 and BAK1, regulating the MAPK activation to induce the plant immune response. In addition, HMGB3 could bind salicylic acid (SA) to inhibit its DAMP activity in Arabidopsis against Botrytis cinerea.[53]

Conclusions

To understand the molecular defense mechanisms that mediate defense of cotton roots against V. dahliae infection, a comparative phosphoproteomic analysis was performed. A total of 639 phosphoproteins, showing significant differences between the mock and inoculated samples in cotton roots in the resistant and susceptible cotton lines, revealed the abundance of known and novel phosphorylated defense-related proteins, signaling-related proteins, transport-related proteins, ATP-binding-related proteins, microtubule-associated proteins, etc. The roles of many of these phosphoproteins are not clearly understood in plant defense. Therefore, this global overview of changes in phosphorylation from cotton roots in response to V. dahliae provides a new perspective for future studies that aim to understand cotton responses to pathogen infection.

Materials and Methods

V. dahliae Culture and Cotton Growth

A defoliating V. dahliae strain Vd080 was isolated from a cotton field in Xinji, Hebei province, China. Then, this strain was grown on PDA medium (potato dextrose agar). The mycelia were collected and cultured in liquid Czapek’s medium as described previously.[54] The final concentration of the spore suspension was adjusted to 1 × 107 conidia/mL with sterile water.

The upland cotton (G. hirsutum) cultivar, Zhongzhimian 2, is highly resistant to V. dahliae, and cultivar Jimian 11 is susceptible to this pathogen. The introgression lines BC6F3:6 were generated from the advanced backcross and repeated selfing using Zhongzhimian 2 as the recipient parent and Jimian 11 as the donor parent. Two materials NIL-R and NIL-S selected from the BC6F3:6 population were resistant and susceptible to verticillium wilt, respectively. NIL-R and NIL-S seeds were delinted with H2SO4 (98%), sterilized in 1% NaClO for 1 min, and then washed three times in sterile water. These seeds were sown in sterile nutrient soil. When the cotyledons were fully expanded, cotton seedlings were grown hydroponically as described previously[55] and maintained in a controlled environment chamber under conditions of 28 °C, 16/8 h photoperiod, and 80% relative humidity.

Inoculation of Cotton Seedlings and Sample Collection

When the first euphylla was fully expanded, the cotton seedlings were inoculated with V. dahliae conidial suspension (1 × 107 conidia/mL) by the dipping root method for 45 min. Control plants were inoculated with liquid Czapek’s medium. Roots were harvested at 1 and 3 days after inoculation (dpi), and the samples were immediately used for extraction of proteins. For gene expression analysis, roots were harvested at 6, 12, 24, 36, 48, 60, and 72 hpi (hours past inoculation) and stored at −80 °C for subsequent use. Three biological replicates were performed for each combination (inoculation/control and time).

Protein Extraction and Tryptic Digestion

Cotton roots were cut into 1 cm segments in length and washed with sterile water four times, each lasting for 5 min. Then, the segments were blotted dry with a piece of filter paper, placed into a mortar, immediately frozen by liquid nitrogen, and grinded. The total protein of roots was extracted by the trichloroacetic acid (TCA)/acetone precipitation method.[56] The Bradford method was used to determine the protein concentration[57] with BSA as a standard. The final protein pellet was resuspended in a lysis buffer and then centrifuged at 20 000g for 30 min. The mixture was reduced by 20 mM DTT for 2.5 h at 37 °C and then alkylated by 40 mM IAA in the dark at 25 °C for 30 min. Then, the final urea concentration was adjusted to 1.5 M by deionized water. The digestion was performed using trypsin with a ratio of 1/50 (w/w) at 37 °C for 18 h. The resulting mixture was desalted by an SPE-C18 column (Waters, WAT051910). The filtrate was collected and freeze-dried. The resulting peptides were stored at −20 °C for further use.

Phosphopeptide Enrichment

The peptide mixture was redissolved in 1.4 mL of precooled IAP buffer (10 mM Na2HPO4, 50 mM MOPS, 50 mM NaCl, pH 7.0) and lyophilized to dryness. TiO2 beads were then added to the peptide mixture with 1X DHB buffer [16% CAN, 0.6% DHB, 0.02% TFA] and then vibrated for 40 min. Then, the supernatant was removed after centrifugation; the beads were put into a tip (Eppendorf, Germany) with a stopper and then washed three times with equilibration buffer I (200 mM NaCl 50% ACN, and 6% TFA) and then three times with equilibration buffer II (0.1% TFA, 30% ACN). Finally, the phosphopeptides were eluted with 50 μL of the elution buffer (NH4OH, pH 10.5) and concentrated under vacuum and then dissolved in 20 μL of 0.1% formic acid (FA) for MS analysis.

Mass Spectrometry Analysis

The mass spectrometry (MS) data were collected by Shanghai Applied Protein Technology, which was performed on a Q-Exactive mass spectrometer (Thermo Scientific) coupled to an EASY-nLC1000 nano LC. The peptides were loaded onto a reverse-phase trap column (RP-C18, 150 μm × 20 mm, Thermo EASY column SC001 traps) connected to the C18 reversed-phase analytical column (RP-C18, 150 μm × 100 mm, Thermo EASY column SC200). Mobile phase A was 0.1% FA (v/v) in acetonitrile–H2O (acetonitrile, 2%), and mobile phase B was 0.1% FA in acetonitrile–H2O (acetonitrile, 84%). A linear gradient was performed at a flow rate of 300 nL/min: 0–220 min, 0–55% B; 220–228 min, 55–100% B; 228–240 min, 100% B.

The mass spectrometer was operated in the positive ion mode with the peptide recognition mode enabled.[29] A data-dependent top10 method was used to acquire the MS data by survey scans (300–1800 m/z) for HCD fragmentation. The resolution for scans was set to 70 000 at m/z 200 and the resolution for HCD spectra was set to 17,500 at m/z 200, and the isolation window was 2 m/z. In each scan cycle, the number of precursors selected for tandem MS was 20, and the charge state screening parameter was set at 2. The underfill ratio was defined as 0.1%, and the normalized collision energy was 30 eV.

Database Search and Quantification

The observed MS/MS spectra were searched against a G. hirsutum L. database (NBI_Gossypium_hirsutum_v1.1.pep.fasta) using the MaxQuant software (version 1.5.3.17)[58,59] identification and quantitation analysis. Trypsin was selected for the enzyme setting. An initial mass tolerance of precursor ions was set at 6 and 20 ppm, and two missed cleavages were allowed. Oxidation of methionine residues was selected as a variable modification, and carbamidomethylation of cysteine residues was selected as a fixed modification. The false discovery rates (FDRs) of peptides and proteins were both set at 1%.

The label-free quantification was performed by MaxQuant software with intensity-based absolute quantification (iBAQ).[58,60] The “Matching Between Runs” feature was enabled and set to 2 min. Only proteins or phosphosites retained were quantified, in at least two of the three biological replicates in at least one sample. The phosphorylated peptides with quantification ratios of >2.0 or <0.5 and p values of <0.05 were defined as statistically different.

Motif-X Analysis

The motifs of identified phosphosites were analyzed using a BioPerl script to extract significantly enriched amino acid motifs.[61,62] There were 13 amino acids in the sequence window, including the modified site and six upstream/downstream amino acids from the modified site. These sequences were used in Motif-x to predict motifs in this study (parameters: width: 13; occurrences: 20; background: species).

Gene Ontology (GO) and KEGG Pathway Analysis

Identified proteins were annotated using a Blast2GO suite.[63] GO annotation (Gene Ontology), GO-Enzyme-Code, and other search tools were used to conduct the phosphoprotein functional classification. The identified phosphoproteins were analyzed against a reference database from P3DB to search for homologues. For annotation augmentation, protein sequences were blasted by InterProScan[64] against EBI (The European Bioinformatics Institute) databases to obtain phosphoprotein domains (motif). Then, GO terms were modulated using the annotation augmentation tool ANNEX.

The KOBAS 2.0 (KEGG Orthology-Based Annotation System) program[65] was applied to map the significant differentially expressed phosphorylated proteins to the terms in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.[66] Due to p values of <0.05, the KEGG pathways were considered significantly enriched in cotton root proteins.

Cluster Analysis of Phosphorylated Peptides

The hierarchical clustering analysis was performed by Cluster 3.0 and Java Treeview software by a Euclidean distance algorithm and an average linkage clustering algorithm for similarity measures and clustering, respectively.

Gene Expression Analysis

Total RNA from roots at 12, 24, 48, and 72 hpi was extracted using a modified cetyltrimethylammonium bromide (CTAB) method.[67] The cDNA was synthesized using the M-MLV reverse transcriptase (TaKaRa) following the manufacturer’s instructions.

The relative expression level of genes was assessed by reverse transcription quantitative PCR (RT-qPCR) based on a previously published method[68] using gene-specific primers (Table S7). The cotton UBQ7 gene was used in RT-qPCR as an internal control. Three independent biological and technical repeats were performed.

Statistical Analysis

All relevant data in this study were presented as mean ± standard error of at least three independent biological replicates. A t-test was carried out to test the statistical significances by SigmaPlot version 12.0 (Systat Software, San Jose, CA).