The seeds of the progenitor soybean cultivar T322 and rn1 mutant (T328H) were planted on germinating papers and incubated under the following growing conditions: 24–26°C day/19–21°C night, 16 h light/8 h dark, and 350 μmol photons m^–2 s^–1 light intensity. Germinating seedlings were phenotyped starting 5 days following sowing for visible necrotic phenotypes. The comparable amounts of necrotic root tissues, as soon as necrotic symptoms started to appear, from the rn1 mutant and healthy tissues from the wild-type T322 plants were harvested and frozen in liquid N₂. The root tissues of random plants were pooled to obtain three replications for each of the wild-type T322 and rn1 mutant plants for proteomics and phosphoproteomics studies.

The proteins from roots of soybean were extracted following the method of Toorchi et al. (2009). Phenol extraction buffer, a mixture of equal volume of extraction buffer and buffered phenol, was prepared just before use. The extraction buffer contains 100 mM Tris–HCl, pH 8.8, 10 mM w/v EDTA, 900 mM w/v sucrose, and 0.4% v/v 2-mercaptoethanol. Buffered phenol was prepared by dissolving phenol in Tris buffer and adjusting the pH to 8.8 with HCl. About 5 g tissue sample was ground to powder under liquid N₂ using a pestle and mortar. Approximately 10 ml phenol extraction buffer was added to the powdered 5 g tissue in a 50 ml tube and protein was extracted by agitating for 30 min at 4°C. The content was centrifuged at 2,800 g for 30 min at 4°C. The upper phase was then transferred to a fresh 50 ml tube. Protease and phosphatase inhibitors were added to 1% concentration. Ice-cold ammonium acetate/methanol solution (100 mM w/v ammonium acetate in 100% methanol) was added to the upper phase in a ratio of 5:1, mixed and incubated at -20°C to precipitate the phenol extracted proteins. The precipitate was collected after centrifugation at 2,800 g for 30 min at 4°C. The pellet was washed twice with of ice-cold ammonium acetate/methanol solution, then twice with ice-cold 80% acetone solution, and finally with ice-cold 70% ethanol. The pellet was solubilized and denatured in a final concentration of 100 mM phosphate buffer pH 7.8, containing 6M urea and 2% CHAPS, 0.2% Triton X-100, 0.1% SDS and 10 mM DTT, sonicated for 5 min and votexed for 0.5 h until completely solubilized. The protein concentration for each sample was determined by Bradford assay, and further quantified by running on a precast NOVEX 12% Tris/Glycine mini-gel (Invitrogen, Carlsbad, CA, United States) along with a series of amounts of E. coli lysates (2, 5, 10, and 20 μg/lane). The SDS gel was visualized with colloidal Coomassie blue stain (Invitrogen), imaged by Typhoon 9400 scanner followed by ImageQuant TL 8.1 (GE Healthcare).

The proteomics and phosphoproteomics analyses were conducted at the Proteomics and Metabolomics Facility, Cornell University to identify the differentially accumulated proteins (DAPs) and phosphopeptides/sites. For TMT labeling, a total of 200 μg protein of each sample was reduced with 9.5 mM tris(2-carboxyethyl)phosphine for 1 h at room temperature, alkylated with 17 mM iodoacetamide for 1 h in the dark and then quenched by additional of 20 mM Dithiothreitol (DTT). The alkylated proteins were precipitated by adding 6 volumes of ice-cold acetone and incubating at -20°C overnight, and reconstituted in 90 μL of 100 mM triethyl-ammonium bicarbonate. Each sample was digested with 18 μg trypsin for 18 h at 37°C. The TMT 6-plex labels (dried powder) were reconstituted with 45 μL of anhydrous ACN prior to labeling and added with 1: 2 ratio to each of the tryptic digest samples for labeling over 1 h at room temperature. The peptides from the following six samples (mut1, mut2, mut3, WT1, WT2, and WT3) were mixed each TMT tag with 126-tag, 127-tag, 128-tag, 129-tag, 130-tag, and 131-tag, respectively. After checking label incorporation using Orbitrap Elite (Thermo-Fisher Scientific, San Jose, CA, United States) by mixing 1 μL aliquots from each sample and desalting with strong cation-exchange (SCX) ziptip (Millipore, Billerica, MA, United States), the six samples in each set were pooled, evaporated to dryness and subjected to a PolyLC SCX cartridge (PolyLC Inc. Columbia, MD, United States) and Sep-Pak C18 cleanup (Waters, Milford, MA, United States) according to the manufacturer’s instructions. The clean samples were then divided to three aliquots: one aliquot with 200 μg peptides used for subsequent high-pH reverse-phase fractionation into 10 fractions as described previously (Yang et al., 2011). The rest two aliquots (500 μg peptides each) were used for enrichment of phosphopeptides.

Enrichment of phosphopeptides by TiO₂ beads: TiO₂ enrichment was conducted using a TiO₂ Mag Sepharose kit (from GE Healthcare) for each 500 μg aliquot. The TMT 6-plex tagged tryptic peptides were reconstituted in 400 μL of binding buffer (1M glycolic acid in 80% acetonitrile, 5% TFA). The TiO₂ slurry (75 μL) was used and incubated with the sample for 30 min at 1,800 rpm vortex. After washing the beads with washing buffer (80% acetonitrile, 1%TFA), the phosphopeptides were eluted with 100 μL of elution buffer (5% ammonium hydroxide) twice. The eluted fraction was dried and reconstituted in 25 μL of 0.5% formic acid for subsequent nano scale LC-MS/MS analysis. The second aliquot (500 μg peptides) was used for sequential enrichment of phosphopeptides by IMAC and TiO₂. The IMAC enrichment was carried out using a PHOS-Select Iron Affinity Gel (Sigma P9740) kit following the vendor-recommended procedure. The sequential enrichment steps were conducted as previously reported (Gao et al., 2017).

The nanoLC-MS/MS analysis was carried out using an Orbitrap Elite (Thermo-Fisher Scientific, San Jose, CA, United States) mass spectrometer equipped with nano ion source using high energy collision dissociation (HCD). The Orbitrap is equipped with a “CorConneX” nano ion source (CorSolutions LLC, Ithaca, NY, United States) coupled with the UltiMate3000 RSLCnano (Thermo, Sunnyvale, CA, United States). The nanoLC-MS/MS was operated in top-15 data-dependent acquisition HCD-MS/MS mode with multiple charged ions above a threshold ion count of 8,000 with normalized collision energy of 37% as described previously for 10 fractions of global proteome (Wang et al., 2014) and two enriched fractions of phosphoproteomics (Gao et al., 2017).

MS raw data files were converted into MGF files for identification and relative quantitation using Proteome Discoverer version (PD) 1.4 (Thermo-Fisher Scientific). A subsequent database search was performed with Mascot Daemon version 2.3 (Matrix Science, Boston, MA, United States) against the Soybean protein RefSeq database (88,647 entries) downloaded on Dec. 1st, 2014, from NCBInr. The default Mascot search settings were as follows: MS peptide tolerance, 15 ppm; fragment mass tolerance, 0.1 Da; trypsin digestion allowing two missed cleavage with fixed modifications of cysteine carbamidomethylation, N-terminal TMT6plex, and lysine TMT6plex; and with deamidation (NQ), oxidation (M), and phosphorylation (STY) as variable modifications. To estimate the false discovery rate (FDR) for a measure of identification certainty in database search, we employed the target–decoy strategy of Elias and Gygi (2007). Specifically, an automatic decoy database search was performed in Mascot by choosing the decoy checkbox in which a random sequence of database is generated and tested for raw spectra as well as the real database. To reduce the probability of false peptide identifications, we considered as identified only those peptides with significant scores (≥25) at the 99% confidence interval by a Mascot probability analysis greater than “identity.” This required that each confident protein identification involve at least two unique peptide identifications indicated in Mascot. In addition, to confidently locate phosphorylation sites, the phosphoRS 3.0 node integrated in PD 1.4 workflow was also used to cross-validate the results of Mascot as reported (Gao et al., 2017).

For protein quantification, we used those proteins that were identified in all six TMT channels with at least two unique peptides in a TMT experiment (Krey et al., 2018). Intensities of the reporter ions from TMT 6-tags upon fragmentation were used for quantification, and the relative quantitation ratios were normalized to median protein ratio for the 6-plex in the dataset. We used Mascot database search engine for processing TMT-based quantitative global proteomics and phosphoproteomics datasets, in which we employed the target–decoy strategy as a measure of identification certainty in database search with 1% FDR as a cutoff. Quantitative comparison of the triplicate samples between the two groups was undertaken with Microsoft Excel software using a student t-test. To identify DAPs related to root-necrosis caused by the rn1 mutation, the quantifiable proteins were tabulated and their differential abundance ratio (treated/control) was log₂ transformed. The log₂ fold values, hereafter referred as fold change (FC), were fit to a normal distribution to obtain the standard deviation of the quantified proteome. The DAPs between rn1 mutant and wild-type 322 were listed using a cutoff of FC value between the rn1 and wild-type 322 based on the criteria of the p-values ≤ 0.05 and FC values ≥ 1.5 or ≤0.67.

The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD032226.

The transcripts of genes encoded the DAPs and phosphorylated proteins were obtained using the fragments per kilobase of exon model per million mapped (FPKM) reads according to the file (genes.fpkm_tracking) listed in the directory of Glycine max Wm82.a2.v1¹ of the Phytozome V13 database (Severin et al., 2010; Goodstein et al., 2012). A stringent filtering criterion of FPKM value 1.0 (equivalent to one transcript per cell; Mortazavi et al., 2008) in the root was used to identify the root-specific expression level of the genes encoding DAPs and phosphoproteins.

The DAPs and phosphorylated proteins were classified into different categories of gene ontology (GOs), biological processes and molecular functions by using GOatools².

Molecular functions of the identified DAPs were annotated by Plant MetGenMap Mercator (Joung et al., 2009) and Kyoto Encyclopedia of Genes and Genome (KEGG) databases (Kanehisa et al., 2017), which were separately performed online³. The pathways of the DAPs including 12 phosphorylated DAPs were generated using MapMan software⁴.

The rn1 mutant exhibited a root necrotic phenotype and enhanced tolerance to the oomycete pathogen P. sojae (Kosslak et al., 1997). To better understand the mechanisms regulated by the Rn1 protein, we compared the root proteomes of the wild-type T322 (Rn1) and mutant T328H (rn1) lines. We used a TMT mass spectrometry-based quantitative proteomics approach to identify the proteins that were differentially accumulated between the lines. Six TMTs were used to label three replicates of necrotic root tissues of the rn1 (T328H) mutant and three of healthy root tissues from the wild-type progenitor T322 plants (Figure 1 and Supplementary Figure 1).

We detected 3,160 proteins, of which 2,180 were found in at least two replications of both necrotic and healthy root tissues. Of the 2,180 proteins, the 150 were DAPs between rn1 and T322 roots, regulated presumably by Rn1, which were selected for further investigation (Table 1). Among the 150 DAPs, 118 showed ≥ 1.5-fold increase and 32 showed ≤ 0.67-fold decrease in the necrotic root tissues as compared to that in the healthy root tissues at p ≤ 0.05. Genes encoding 138 of the 150 DAPs showed detectable expression levels (Transcript Number ≥ 1) in roots of the soybean cultivar Williams 82 (Table 1).

To understand the role of the 150 DAPs in the spontaneous cell death process leading to root necrosis, DAPs were classified based on biological processes, cellular components and molecular functions by using the plant GO database (Harris et al., 2004). As shown in Figure 2, the DAPs were grouped into 17 basic biological processes including cellular process, metabolic process, single-organism process, response to stimulus, and biological regulation. According to cellular components, the DAPs were grouped into 12 types including cell part, cell, organelle and membrane. Based on molecular functions, the DAPs were grouped into eight classes including binding, catalytic activity, and transporter activity.

Based on the plant-specific database of Clusters of Orthologous Groups of proteins (COGs; Tatusov et al., 2003), the 150 DAPs were classified into 17 groups (Supplementary Figure 2). As expected, there were more upregulated proteins than down-regulated ones among 12 of the categories. In five categories, viz. (i) nucleotide transport and metabolism, (ii) carbohydrate transport and metabolism, (iii) intracellular trafficking, and vesicular transport, (iv) chromatin structure and dynamics, and (v) translation, ribosomal structure and biogenesis, more down-regulated proteins than upregulated proteins were represented in the rn1 mutant compared to T322 wild-type line suggesting that transcription and translation machineries may be suppressed in the necrotic root tissues.

To understand how the rn1 mutation affects the accumulation of proteins among the known pathways, we analyzed the proteins using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. This database is an ideal tool to predict the gene function by linking genomic information with higher order functional information of cellular processes including metabolism, membrane transport, signal transduction and cell cycle (Kanehisa et al., 2017). The overall enrichment of the 150 DAPs is shown in Figure 3A. The majority of the DAPs are involved in metabolism (46.71%). Accumulation patterns of a large number of uncharacterized soybean proteins (27.63%) and with unknown function (15.79%) are also affected by the rn1 mutation. DAPs were classified to 42 pathways under five categories with known functions (Figure 3B). Several of the proteins are involved in primary and secondary metabolism. Pathways enriched among the DAPs include glyoxylate and dicarboxylate metabolism, propanoate metabolism, cyanoamino acid metabolism, sulfur metabolism, arachidonic acid metabolism, alpha-linolenic acid metabolism, linolenic acid metabolism, carbon fixation in photosynthetic organ, terpenoid backbone biosynthesis, isoflavonoid biosynthesis, and phenylpropanoid biosynthesis. DAPs were mapped to pathways using MapMan bin codes (Supplementary Figure 3). The 21 DAPs were found to be associated with RNA-protein syntheses pathways (Supplementary Figure 4).

The rn1 mutation enhances tolerance to the oomycete pathogen P. sojae in soybean (Kosslak et al., 1997). Plants respond to pathogen attacks by a rapid change in gene expression levels leading to accumulation of biotic stress-related defense proteins and metabolites or compounds. MapMan was used to analyze the 150 DAPs to determine which DAPs are involved in pathways that are activated by biotic and abiotic stresses. Among the 150 DAPs, 53 and eight DAPs were found to be involved in biotic and abiotic stresses, respectively, (Figure 4).

Of the identified biotic stress-related proteins, five DAPs (Glyma.01g126600, Glyma.03g147700, Glyma.08g341300, Glyma.08g341400, and Glyma.16g212400) are pathogenesis-related (PR) proteins are highly induced in the rn1 mutant as compared to that in T322. They are also expressed in roots and could be involved in enhancing tolerance to P. sojae (Figure 4).

Several DAPs responsive to biotic stress are involved in hormone signaling, redox reaction, signal transmission, cell wall modification, and synthesis of secondary metabolites observed in plant–pathogen interaction process. For example, 11 DAPs, auxins-associated proteins (Glyma.15g115600 and Glyma. 06g260200), ethylene-associated proteins (Glyma.02g048400, Glyma.02g268000, Glyma.05g222400, and Glyma.02g048400), and JA-associated proteins (Glyma.13g347600, Glyma. 15g026300, Glyma.13g109800, Glyma.01g235600, and Glyma.11g007600), participate in hormone-regulated defense signaling pathways and mediate the biosynthesis of plant defense compounds (Figure 4).

Protein phosphorylation is a key mechanism involved in regulating protein functions. In this study, we identified a total of 573 phosphopeptides of 146 proteins through LC-MS/MS analysis (Supplementary Table 1). Of the 573 phosphopeptides, only 234 were unique phosphopeptides. Of the 234 peptides, only 165 peptides were identified. GO analysis of the phosphorylated 146 proteins assigned the phosphoproteins (PhosPs) to 19 classes based on biological processes, 11 based on cellular components, and 10 based on molecular functions (Supplementary Figure 5). Based on molecular functions, PhosPs were grouped into many classes of proteins with binding and catalytic activity, transporter activity, nucleic acid binding, signal transducer activity, electron carrier activity, molecular transducer activity, molecular function regulator and antioxidant activity. PhosPs with antioxidant activity could be involved in conferring disease resistance.

Among the 146 PhosPs, 143 proteins displayed detectible expression levels in root (Supplementary Table 1). Only two of 143 proteins were found to be significantly changed in their steady state levels (≤0.67-FC, p-value ≤ 0.05) due to the rn1 mutation. One is Phosphoglucomutase/phosphomannomutase family protein (Glyma.05G237000.1), and the other one is an unknown protein (Glyma.16G037900.1). The levels of phosphorylation for these two differentially accumulated phosphoproteins are significantly reduced in the rn1 mutant root tissues as compared to that in the healthy root tissues of T322.

We analyzed the proteomic and phosphoproteomic data to investigate how the rn1 mutation contributes to the phosphorylation statuses of the PhosPs. For this study, we selected 88 PhosPs that contain phosphopeptides with significantly different levels of phosphorylation between the root tissues of rn1 and T322. For each PhosP, nine comparisons were made for “levels of phosphorylation” (A) and “levels of PhosP” (B); i.e., each of the three replications of the rn1 root sample was compared to each of the three root samples of the T322 wild-type plants for phosphorylation levels (A = rn1/T322) or steady state protein levels of PhosP (B = rn1/T322). The “A” values from nine comparisons were then compared to respective nine “B” values from nine comparisons of PhosP levels to obtain nine A/B ratios for each of the 88 PhosPs.

The phosphoproteomics analysis identified 22 phosphopeptides, the levels of phosphorylation of which were significantly different between rn1 and T322 lines (Table 2). Among the 22 phosphopeptides, two showed enhanced phosphorylation levels, and 20 showed reduced phosphorylation levels in the rn1 compared to the wild-type 322. The PhosPs carrying these 22 phosphopeptides were grouped into 14 classes based on basic biological processes, 11 based on cellular components and seven classes based on molecular functions (Figure 5A). In particular, some of these PhosPs are involved in the immune response process.

The 22 PhosPs were enriched into 12 pathways (Figure 5B). It is clear that some proteins are related to the oxidative phosphorylation and pentose phosphate pathway, viz. (i) H⁺ ATPase (Glyma.09G056300) involved in oxidative phosphorylation (Supplementary Figure 6); (ii) phosphoproteins related to the pentose phosphate pathway enzymes are glucose-6-phosphate dehydrogenase 6 (Glyma.16g063200 and Glyma.19g082300), phosphoglucomutase/phosphomannomutase (Glyma.08g044100 and Glyma.05g237000) family protein, (Supplementary Figure 7). MapMan analysis revealed that seven of the 22 PhosPs are induced by biotic and abiotic stresses (Figure 6). Especially, phosphorylation levels of two highly similar type II metacaspases (Glyma.08G233300 and Glyma.08G233500) were significantly reduced in the rn1 mutant root tissues as compared to that in the healthy root tissues of T322, suggesting that the dephosphorylation or reduced phosphorylation of type II metacaspases may contribute toward initiating spontaneous cell death observed in roots of the lesion mimic rn1 mutant (Supplementary Figure 8). The phosphopeptides of the Glyma.08G233300 and Glyma.08G233500 with significant differences in phosphorylation levels (Table 2 and Supplementary Table 2), and Glyma.15G219100.1 (Table 3 and Supplementary Table 2) with statistically non-significant reduced phosphorylation level were localized to the C-terminal region of the P20 caspase-like domain (Supplementary Figure 9). The MS/MS spectra of the identified phosphopeptides from all three type II metacaspases (Glyma.08G233300, Glyma.08G233500, and Glyma.15G219100) along with their quantitative changes between the rn1 mutant and wild-type 322 lines are shown in Figure 7.

We have identified 24 less abundant PhosPs containing phosphopeptides that exhibited significantly different phosphorylation levels between the rn1 mutant and T322 wild-type root tissues (Table 3 and Supplementary Table 2). The steady state levels of these peptides could not be detected because of their low abundance. Therefore, phosphorylation levels of the 24 phosphopeptides between rn1 mutant and T322 lines could not be compared as in for 22 PhosPs (Table 2), and are considered as candidates for future studies. These phosphopeptides were classified into 15 biological processes, 11 cellular components, and three molecular functions based on GO terms (Supplementary Figure 10A). The 24 phosphopeptides can be classified into nine pathways (Supplementary Figure 10B). MapMan analysis of the 24 phosphopeptides revealed that four proteins (Glyma.08G243000, Glyma.18G195500, Glyma.05G011600, and Glyma.08G297500) involved in biotic and abiotic stresses pathways (Supplementary Figure 11).